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Abstract. A mathematical framework for Continuous Time Finance based on operator al- 
gebraic methods offers a new direct and entirely constructive perspective on the field. It also 
leads to new numerical analysis techniques which can take advantage of the emerging massively 
parallel GPU architectures which are uniquely suited to execute large matrix manipulations. 

This is partly a review paper as it covers and expands on the mathematical framework un- 
derlying a series of more applied articles. In addition, this article also presents a few key new 
theorems that make the treatment self-contained. Stochastic processes with continuous time 
and continuous space variables are defined constructively by establishing new convergence es- 
timates for Markov chains on simplicial sequences. We emphasize high precision computability 
by numerical linear algebra methods as opposed to the ability of arriving to analytically closed 
form expressions in terms of special functions. Path dependent processes adapted to a given 
Markov filtration are associated to an operator algebra. If this algebra is commutative, the 
corresponding process is named Abclian, a concept which provides a far reaching extension of 
the notion of stochastic integral. We recover the classic Cameron-Dyson-Feynman-Girsanov- 
Ito-Kac-Martin theorem as a particular case of a broadly general block-diagonalization algo- 
rithm. This technique has many applications ranging from the problem of pricing cliquets 
to target-redemption-notes and volatility derivatives. Non-Abelian processes are also relevant 
and appear in several important applications to for instance snowballs and soft calls. We show 
that in these cases one can effectively use block-factorization algorithms. Finally, we discuss 
the method of dynamic conditioning that allows one to dynamically correlate over possibly 
even hundreds of processes in a numerically noiseless framework while preserving marginal 
distributions. 
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1. Introduction 



The goal of this paper is to attempt to consolidate and present a number of mathematical 
methods developed over several years by myself and collaborators while addressing concrete 
problems in derivative pricing theory. The results scattered across a number of papers which 
are collected here have been complemented with a rigorous ab initio treatment and a few key 
theorems which make the framework mathematically self-contained. This results in a quite 
comprehensive approach to the theory of Stochastic Processes and Mathematical Finance which 
is novel in that it is fully constructive and perhaps has applications beyond the realm of Financial 
Engineering. 

There are several traditions of Constructive Mathematics. One attempts to re-derive classical 
results of real and functional analysis based on a restrictive constructivist logic according to 
which no mathematical object can be considered unless one specifies explicitly how to construct 
it, see (Bridges and Richman 1987) and (Bishop 1967). Along another tradition, Constructive 
Field Theory, see (Glimm and Jaffe 1987), aimed at establishing the existence of interacting 
quantum field theories by providing a constructive procedure for computing n-point functions and 
demonstrating that they satisfy a set of axioms. Measure theoretic probability and the related 
theory of stochastic processes, (Doob 1953), does not seem to be understandable constructively. 
The PDE approach in (Feller 1961) and the harmonic analysis approach in (Bochner 1955) are 
instead essentially constructive but do not delve into the theory of stochastic integrals and path 
dependent processes and into lattice discretization schemes. 

The main motivation that guided this research is the creation of an engineering framework for 
exotic financial derivatives. Efficient computability on current hardware has been and remains 
throughout this article our key motivating concern. To this end, we work towards an algebraiza- 
tion of Probability Theory that reduces all calculations to matrix manipulations which can 
be performed efficiently and in particular to matrix multiplications. Similarly to the standard 
framework of algebraic topology, (Spanier 1966), we consider processes taking values in separable 
topological spaces and approximate continuous domains by means of simplicial sequences. To 
establish convergence in the continuous limit, we directly estimate convergence rates for proba- 
bility transition kernels in the continuous space limit following an approach similar in spirit to 
Constructive Lattice Field Theory, see (Glimm and Jaffe 1987). Similarly to constructive field 
theory, sets of axioms on n-point functions are used to identify processes and renormalization 
group transformations are used to control the continuous limit. Following (Naimark 1959), the 
approach is grounded upon the algebraic theory of integration on locally compact Hausdorff 
separable topological spaces. 

Calculations with stochastic processes are carried out using operator methods developed in 
Quantum Mechanics, (Landau and Lifshits 1977) and systematized in Mathematical Physics 
references such as (Reed and Simon 1980). In Finance, operator methods have been developed 
along two independent and non-overlapping streams of research, one by Ait-Sahalia, Hansen 
and Scheinkman who focused on econometric estimations in a series of papers reviewed in (Ait- 
Sahalia et al. 2005), see also (Ait-Sahalia 1996), (Hansen et al. 1998), (Hansen and Scheinkman 
1995). The second stream of research is by the author and collaborators who instead worked 
on derivative pricing for path dependent and correlation derivatives, see (Albanese et al. 2005- 
2006&), (Albanese and Chen 2004a), (Albanese and Chen 2004&), (Albanese and Kusnetsov 2005), 
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(Albanese and Lawi 2004), (Albanese et al. 2006a), (Albanese and Trovato 2006), (Albanese 
and Trovato 2005), (Albanese and Vidler 2007), (Albanese and Jones 2007) and (Albanese and 
Osseiran 2007). In this paper, we attempt to systematize the mathematical framework of pricing 
theory in the operator formalism from our own viewpoint, reserving to future work the task of 
pursuing overlaps with the econometric literature. 

The references quoted above are all relevant to our undertaking and provided motivations 
on many levels. However, in an effort to keep this writing self-contained, we are not going to 
assume any previous knowledge of the reader. 

To ground the mathematical framework, we obtain sharp pointwise convergence estimates 
for probability kernels and its derivatives. More precisely, we show that probability kernels 
converge pointwise at rates of order 0(h 2 ), where h is the lattice spacing. The result applies to a 
large class of diffusion processes and extensions thereof including smooth space inhomogeneities, 
regime switching, finite activity jumps and some degree of time inhomogeneities. We also prove 
similar convergence results for the fast exponentiation method, our preferred numerical method 
for exponentiating Markov generators, showing that errors in this scheme are also of order 
0(h 2 ), in the sense of pointwise convergence for the probability kernel. In the particular case of 
Brownian motion, we show that similar 0(h 2 ) pointwise error estimates apply also to derivatives 
of the probability kernels and that the power 2 in the 0(h 2 ) bounds is actually sharp. 

The interest in convergence estimates was prompted by the desire of understanding the mech- 
anisms behind the empirically observed smoothness and robustness in the calculation of price 
sensitivities with the methods in our applied papers. We also observed empirically that the 
fast exponentiation algorithm is stable under single precision floating point arithmetics. We find 
that key to a high precision numerical framework for sensitivities is handling the time coordinate 
either as continuous or very finely discretized. A sufficiently fine discretization is defined as one 
for which explicit differentiation schemes are stable and typically correspond to a hourly time 
scale in applications. Typically, weekly time steps would permit stable implicit differentiation 
schemes but would not allow for as much stability in the calculation of price sensitivities and the 
probability kernels we require to evaluate and manipulate. This motivates us to avoid implicit 
differentiation schemes on coarse time intervals. 

Measure changes and time changes are defined constructively and a version of the Fundamen- 
tal Theorem of Finance is re-obtained. The Cameron-Martin-Girsanov's theorem, see (Cameron 
and Martin 1949), and Ito's lemma, see (Ito 1949), are proved twice in different ways with oper- 
ator methods. We also derive the Feynman-Kac formula, see (Feynman 1948) and (Kac 1948). 
One of the key results is an extension of the Feynman-Kac-Ito formula in three different direc- 
tions. This formula concerns the characteristic function of a stochastic integral over a diffusion 
process. In our formalism, this formula becomes a block-diagonalization algorithm for large 
matrices associated to path-dependent processes. The extension we discuss (i) covers Markov 
processes more general than diffusions, (ii) allows for a class of path-dependent processes we 
name Abelian which extends the notion of stochastic integral and (iii) generalizes the harmonic 
analysis framework to include extensions of trigonometric Fourier transforms. The theory of 
Abelian processes finds numerous practical applications to path dependent options and is appli- 
cable to the great majority of path-dependent payoffs, from volatility swaps to cliquets, range 
accruals, lookback options, target redemption notes and more. We also give a version of Dyson's 
formula to accelerate the pricing of path-dependent options given by Abelian processes by means 
of a moment expansion. Non-Abelian processes are more difficult to handle but we single out 
a class admitting block-factorizations (as opposed to a block-diagonalizing transformation) and 
which are also amenable to numerical analysis by matrix algebra. Finally, we illustrate the 
method of dynamic conditioning that allows one to correlate possibly numerous processes by 
means of kernel manipulations while preserving marginals and not incurring into dimensional 
explosion. 

The mathematical methods in this article are particularly efficient as they lend themselves to 
transparent hardware acceleration on the emerging multi-core GPU hardware platforms. These 
massively parallel architectures are based on low-cost technologies that have been developed 
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for the games and high definition markets and arc uniquely suited to implement BLAS Level 3 
routines such as matrix-matrix multiplications with high efficiency. See also (Goto and van de 
Geijn to appear) for a state of the art account on matrix-matrix multiplication software on 
CPUs. 

The paper is organized as follows. In Section 2, we introduce the notion of simplicial sequence 
which is key to devising approximation schemes for continuous valued process by means of a 
sequence of Markov chains. In Section 3, we consider a general definition of path functional and 
in Section 4 we give a general description of a stochastic process by means of an n-point function. 
Markov Processes are introduced in Section 5, martingales and monotonic processes in Section 6. 
In Section 7, we derive the Fundamental Theorem of Arbitrage Free Pricing Theory. The classical 
results on weak convergence of Markov generators by Bernstein, Bochncr, Kyntchine and Levy 
are re-obtained in Section 8. Time homogeneous Markov Processes, fast exponentiation and 
spectral methods are described in Section 9. In Section 10, we carry out a constructive analysis 
of Brownian motion and prove convergence estimates. In Section 11, we study the spectrum of 
diffusion generators. Sharp pointwise kernel estimates are extended to general diffusion processes 
in Section 12. In Section 13, we give estimates for the convergence rate of time discretisation 
schemes of the type we advocate for applications, i.e. based on fast exponentiation. Section 
14 reviews the derivation of hypergeometric Brownian motion and particular cases such as the 
CEV model. In Section 15, we study stochastic integrals and obtain Ito's formula for diffusion 
processes and of Girsanov's theorem. Section 16 contains a derivation of the Feynman-Kac 
formula for bridges over general Markov processes. The general notion of Abelian process in 
continuous time is introduced in Section 17. In Section 18, we discuss the discrete time case. In 
this section, we introduce also the notion of non-resonant block-diagonalisation scheme which 
provides a numerically useful extension of Fourier analysis based on trigonometric functions. 
Dyson's formula and moment expansions arc in Section 19, covering the uni-variate case, and 
Section 20, covering the multivariate case. These two sections include applications to exotic 
volatility derivatives. Block factorizations and applications to snowballs and soft calls are in 
Section 21. Dynamic conditioning and multi-factor correlation modeling is discussed in Section 
22. Conclusions end the paper. 

2. Measure Theory on Simplicial Sequences 
Let d > be an integer, consider the space M d and the sequence of lattices h m Z d where 

Definition 1. (Lattices.) If A C h m Z d , the convex hull of A in M d is denoted by Hull(A). 
The interior of A C h m Z d is denoted with Int(.A) and is defined as the set of all sites x £ A 
contained in A along with each one its neighbors in h m Z d at distance h m . 

Definition 2. (Simplicial Sequences.) A bounded simplicial sequence is given by an integer 
too > 0, and a sequence of subsets A m £ h m T, d defined for all to > too such that 

• Hull(^4 m ) C Hull(A m /) whenever to < to' '. 

• for all to > and all internal points x £ lnt(A m ) there is a M > and an e > such 
that, for all m' > M and for all y £ h m /Z d with d(x, y) < e we have that y £ A m > . 

We follow a constructivist logic paradigm according to which in order to identify a set or 
a sequence of sets one has to explicitly state how to construct it, possibly with a recursive 
algorithm, and it must be possible to decide each step of the recursion in a finite number of 
logical steps. 

Definition 3. (Lattice Functions.) Let A = (A m ), m > too be a bounded simplicial sequence. 
A real valued simplicial function, denoted by f : A — > M. is defined as a sequence of functions 
fm ■ A m — > R such that f m '{x) = fm{x) f or all m' > to > m and all x £ lnt(A m ). The function 
f : A — > R is said uniformly bounded if there is a constant c > such that |/ m (a;)|< c for all 
x £ A m . The function f : A — ► M is uniformly continuous if it is uniformly bounded and for all 
e > there is a 5 > such that if x,y £ h m Z and d(x, y) < S we have that \f(x) — f(y)\< e. 
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Definition 4. (Equivalence of Lattice Functions.) Let A = (A m ), m > mo, be a bounded 
simplicial sequence and f — (/,„) and g = (g m ) are two uniformly continuous real valued func- 
tions on A. We say that these series provide equivalent representations of the same function in 
case, for all m > and all all x G h m 7L we have that lim m /^ 00 |/ m /(x) — g m >(x)\— 0. 

Let C(A) be the set of all continuous functions on A endowed with the natural structure of C* 
algebra given by the operations of sum, multiplication by a scalar and pointwise multiplication. 

Definition 5. (Integrals.) Let A = (A m ), m > mo, be a bounded simplicial sequence. An 
integral is given by a sequence I = (I m ) where I m is a linear functional on the linear space of 
functions f. m : A m — ► K such that 

• Im{fm) > whenever f m (x) > for all x e A m . 

• The following limit exists for all continuous functions f — (f m ) e C(A): 

(2.1) /(/) = lim I m (f m ) 

m—>oo 

• The junctionals I m (fm) defined above satisfies the following bound: 

(2.2) I m (f m ) < c||/ ro |U 

for some constant c > and all functions f m . 
An integral is said to correspond to a probability measure if I m (l) = 1. 

Definition 6. (Equivalent Integrals.) Let A = (A m ), m > mo, be a bounded simplicial 
sequence and let (I m ) and (J m ) be two integrals. We say that these series provide equivalent rep- 
resentations of the same integral I i/lim m _ >00 |/ TO (/ TO ) — J m (/ TO )|= for all uniformly continuous 
functions (f m ). 

Given an integral I on the simplicial sequence A, for all p > 1, one defines the following 
semi-norm on the space of continuous functions C(A): 

(2-3) \\f\\ P =(I(\f\ P )) 1/p - 

Let L P (A;L) be the linear space obtained by completing C(A) with respect to the semi-norm 
1 1/| | j,p and identifying equivalence classes of functions at zero distance. More precisely, L P (A; L) 
is the linear space of the Cauchy sequences (f m ) in C(A) with respect to the norm ||/||/, p modulo 
the linear space of the Cauchy sequences converging to a limit of zero L p -norm. 

Definition 7. (Summable Functions.) A function is called summable if it is in L 1 (A;L) and 
square summable if it is in L 2 (A; I). 

Theorem 1. (Monotonic Sequences.) Let fk be a non-decreasing sequence of summable 
functions and consider the function f = lirrifc^oo fk and the limit I = ]im.i c — >00 I (fk) ■ Then 
either I = oo or f is summable and I = 1(f). 

Proof. If I < oo, then whenever k > j we have that ||/fe — /j||i= I(fk — fj)- Since the sequence 
I(fk) is uniformly increasing, it is Cauchy. Hence the sequence fk is Cauchy in L 1 (A;I). Its 
limit / therefore belongs to L 1 (A; I) since this space is complete. □ 

Definition 8. (Measurable Functions.) A function is called measurable if it can be repre- 
sented as the limit of a non- decreasing sequence of summable functions. 

Theorem 2. (Dominated Convergence.) Let fk be a sequence of measurable functions, 
suppose there is a summable function g such that |/fe(a;)|< \g(x)\ and suppose also that the 
limit f(x) = linifc^oo fk(x) exists in a pointwise sense. Then f is a summable function and 
1(f) = lim fc _ >00 /(/ fe ). 

Proof. Consider the sequence Uk constructed iteratively so that Ui = min(ui,/) and Uk — 
max(ufc_i, min(/fe, /)). The sequence Uk is uniformly non-decreasing and converges to f(x). 
Furthermore, /(|tifc|) < I(\g\) for all k. Hence / is summable and /(/) = limfe^oo I(fk). □ 



6 



CLAUDIO ALBANESE 



Let V be the linear space of all functions which can be expressed as the limit of a non- 
decreasing sequence of continuous functions f(x) — lim^oo fk(x) such that the following norm 
is finite 

(2.4) ||/|U= lim H/felU. 

k — >oo 

L°° (A; I) is defined as the completion of the linear space V with respect to the uniform norm. 

Definition 9. (Essentially Bounded Functions.) A function is called essentially bounded if 
it is in L°°(A;I). 

Definition 10. (Absolute Continuity.) Let I and J be two integrals on the simplicial sequence 
A = (A m ). I is said absolutely continuous with respect to J if these two integrals admit simplicial 
representations I = (I m ) and J = (J m ) and there exists a J—summable function g = (g m ) such 
that I m (fm) = J m (g m f m ) for all I-summable functions f = (f m ). 

3. Path Functionals 

To introduce processes we need to specify the notion of measure on a set of paths. One 
could possibly introduce simplicial sequences in space-time but we refrain from doing so and 
consider instead the time coordinate as a continuous variable. In the time direction, we are only 
going to be using Ricmann integrals for piecewise smooth functions, so that the underlying time 
discretisation one might possibly imagine is quite straightforward and keeping track of it would 
only lead to notational complexities. 

Let us consider a bounded simplicial sequence A — (A m ) and let [T, T'] C R be a fixed finite 
time interval. Let us consider the sequence A = (A m ) where A m = A m x [T,T']. 

Given an integer q > 0, let denote the set of all functions y. : [T,T'] — > A m which are 
constant on a family of q mutually disjoint sub-intervals of the form [t/j 4 , tfc i+1 ), i = 0, ..q with 
t ka = and t kq+1 — T' , spanning the interval [T,T']. Notice that (H q n ) for q fixed is itself a 
finite dimensional manifold with boundaries. In fact, a function y G is characterized by the 
ordered sequence to — T < t\ < .... < t q+ i = T' of time points between which y t is constant and 
a set of values y , such that, y s = yi for all s e [t i} t i+ i) and for all i = 0, ..q — 1. 

The path spaces for q = 0, 1, ... can be regarded as nested into each other Ti.o(T,T') C 
Wi(T,T') C .... In fact, a path y t in H q m is also a path in H m >{T,T') if m < to'. Let H(T,T') = 
U m =o4,..'^m be the union of all path spaces containing paths with a finite number of jumps. 

Definition 11. (Function Algebras.) Let us denote with C q (A), the C* algebra of all uni- 
formly continuous sequences of simplicial functions F = F m : 7i q n —> K endowed with the 
operations of sum, multiplication and with respect to the uniform norm defined as follows: 

(3.1) ll*1loo= I™ sup \F(y.)\. 

Definition 12. (Path Functionals.) A path functional is a sequence F = (F q ), q = 0,1, ... of 
continuous simplicial functions F q — (FfiA <E C q (A) satisfying the following mutual compatibility 
condition: 

(3.2) F£{y.) = Fl{y.), for all y.eH q n and all q' > q. 

Definition 13. (Non-anticipatory Path Functionals.) Let F(y.,t) = (F q (y.,t)) t e [T,T'\ 
be a one-parameter family of path functionals. One says that this is a non anticipatory path 
functional if 

(3-3) F q m (y.,t) = F q m (y[,t) 

whenever y s = y' s for all s < t. 

Intuitively, non-anticipatory path functionals are indifferent to information about the re- 
alization of the path y. in the argument at future times s > t. An elementary example of 
non-anticipatory path functional is given by a function of two arguments 

(3-4) Fl(y.,t) = F (y t ,t) 
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where Fq : Aoo x [T, T'\ — > K is a piecewise smooth one-parameter family of functions. Ap- 
plications typically call for functions F Q (jj 7 t) which are piecewise smooth in t for each y G A 
with possibly a discrete set of jump discontinuities. We follow the usual convention according 
to which if jumps occur, then the discontinuity is of cadlag type, i.e. right continuous and with 
a left limit. 

A less elementary example of non- anticipatory path functional one often encounters is given 
by integrals of the form: 

(3.5) F^(y.,t)= / ds 1 ... ds n a(t;s 1 ,....s k )F 1 (y Sl ,s 1 )...F k (y Sk ,s k ). 

Jt Js k _ 1 

where the Fj : x [T, T'] — > R, j = 1, ...k, are one-parameter families of lattice functions and 
also a(t; So, ....Sfc) is a function of the time coordinates. Applications typically call for functions 
Fj(y,t) which are piecewise smooth in t for each y E A with possibly a discrete set of jump 
discontinuities. Also in this case, we follow the convention according to which at the points of 
jump discontinuity the function is cadlag. The same regularity assumptions will be postulated 
for the functions a(t; sq, ....Sk) with respect to each of its arguments. Although the regularity 
assumption for these path functionals is sufficient for applications, they are not strictly needed 
and can be relaxed by taking limits. 



4. n-point Functions 

Definition 14. (Filtered Probability Spaces.) Consider the C* algebra V generated by 
the functionals of form (3.4-) and (3.5) by taking linear combinations of finite products and 
completing the resulting normed space. A filtered probability space upon the lattice A is defined 
as a bounded linear functional on the C* — algebra V . 

A constructive definition of a stochastic process can be given with various degrees of gener- 
ality. We don't aim here at the utmost generality but instead at pedagogical simplicity and at 
ultimately explaining how to frame the theory of Markov processes. As a step toward this goal, 
we introduce here a family of time ordered n— point functions 

(4.1) $m(yuti;—;y n >tn\yo) 

defined for yi, ...y n £ A m and t\, ...t n £ [T,T'] which is assumed to be piecewise differentiable 
in the time coordinates. In addition, we assume that the following properties hold for all T < 
h < ... < t n < T: 

(SP1) $(2/i,*i; ...\y n ,t n \y ) > Vy , Vi, -Vn € A, 

(SP2) ®(yi,h;...;y n ,t n \y ) = l Vj/ € A, 

yi,...y n £A 

(SP3) ...;yi,U;y i+1 ,U +1 ; ...;y n ,t n \y ) = 

if ti — for some i = l...n — 1 and if yi ^ yi+i 

(SP4) y i ,t i ;y i+1 ,t i+1 ; ...;y n ,t n \y ) 

v 

$(yi,ti; -;yn, t n \yo) Vj/i, ..y n e A. 

Notice that for each fixed starting point yo and fixed sequence T < t\ < ... < t n < T', the 
function <i> m (yi,ii; ...; y n , t n \yo) is a probability distribution function in the arguments yi, ...y n . 
This is interpreted as the probability distribution density for paths starting from the site y^ at 
time T and achieving the values y\, ...y n at times t\, ...tn. 

Definition 15. (Stochastic Processes.) An adapted (stochastic) process is given by a non- 
anticipatory path functional F € V and a measure on P defined by a family of n— point functions 
(•!•)„ 1.2... 
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Given an adapted process of the form F t = F(y t ,t) E V and given a t > T, the expectation 
subject to the initial condition yx — ya of the value attained by the process F t at time t is given 
by 

(4.2) E[F t \y T = y ]=J2 $(yx,t\y )F(y u t). 

yi eA 



For an adapted process of the form (3.5 1 instead, the conditional expectation is given by 
(4.3) 

E[F t \yT = yo] = Y / ds 1 ... ds n ^(yi,si]...;y n ,s n \y )a(t;si,...s n )Fi(yi,si)...F n (y n ,s n ). 
//:..., Jr J - 

Also higher moments can be computed. To keep expressions simple, consider a path functional 
of the form 

(4.4) F t = f ds 1 a(t;s 1 )F 1 {y Sl ,s 1 ). 



The variance of this process at a future point in time conditional to the starting point at time 

T is given by 

(4.5) 

E \ F t\VT = Vo] = 2 ^ / dsi ds2$(yi,s 1 ;y2,S2\yo)a(t;s 1 )a(t;s 2 )F 1 (y 1 ,s 1 )F 1 (y2,S2)- 



Notice that the factor 2 compensates for the time-ordering needed to recast the expression in 
such a way that we can apply to it a 2-point function. This expression and its generalizations 
are used extensively in the applications to path-dependent options discussed below. 

More generally, consider two path functionals F t and G t G L 1 (W(T, T'),fi), where F t is defined 



as in (3.5 1 and 



(4.6) G t = / dsi... dSqb(t;s 1 ,....s q )Gi(y Sl ,s- l _)...G q (y Sg ,s q ). 



Then we can compute a mixed moment as follows: 



E[F t G t \vT =Vo] = y~] / dsx... I ds n+q $(yi, sr, y n + q , s n+q ) 

^7r(l) ) ^7r(n) 

)^i(y77(i),s w (i ) )...i ;i „(2/ 7r („), s l(n) ) 



E 



&(£; s^^j.).!), ••■•s x („ +g ))Gi(?/ 7T -( n+1 ), s T („+i))...G 9 (y 7r ( n+I j), S7r(,i+ 9 )), 

where the sum ranges over all the permutations ir of the time ordered sequence Si < ... < s n+q 
such that 5^(1) < ... < s^n.) and s 7T („+i) < ... < S7r(n+<j). Higher conditional moments of any 
path functional of the form above can be evaluated in a similar way. 

n-point functions Si\ y n , s n \yo) are conditioned to the starting point j/o- It is straight- 
forward to define variations of these n-point functions which reflect conditioning on an initial 
stretch of a path y s for s £ [T, t) , where t > T. Conditioning to a past history is equivalent 
to restricting the integral over each manifold Ti(T, T') to a sub-manifold which is part of its 
boundary. In general, this results in rather clumsy expressions which are difficult to compute. 
In the next section we specialize to the Markovian case for the underlying lattice process and in 
this case there are no memory effects and conditioning is more straightforward. 

Definition 16. (Radon-Nykodim Derivative.) Consider two sequences of n-point functions 
on the simplicial sequence A m $^(yi, h; y„, t„\y Q ) and $^(2/1, h; y n , t n \y ) , where m = 
1,2, .... The Radon-Nykodim derivative of < & 1 with respect to <i> 2 is given by the path functional 
defined as follows: 

(4.7) ,(„.)= Hm 

n-*oo $4(y tl ,t i; ...;y tn ,t n \y ) 
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where U = T + ^ (T" — T) . Notice that the Radon- Nykodim derivative might possibly be infinite 
on some paths. 

The measure in path space given by is said to be absolutely continuous with respect to 

($>rn) if th- e Radon-Nykodim derivative is finite and summable. 

Finally, a word on stopping times. 

Definition 17. (Stopping Times.) A stopping time is an adapted process r t = r{y.,i) which 
can take only two values, by convention and 1. 

Stopping times are often used in conjunction with other adapted processes F t to construct 
stopped versions of it. If F t — F(y.,t) is an adapted process and Rt — R(y.,t) is a second 
process, then a stopped version of F t corresponds to the adapted process of function F(y.,t), 
where F(y. , t) = F(y. , t) if r t = and F(y. , t) = R{y. , t) if r t = 1. 

5. Markov Processes 

Definition 18. (Markov Propagator.) A Markov propagator on the simplicial sequence 
A m ,m > mo is defined as a sequence of functions U m {yi,ti;y 2 ,t 2 ) where 2/1,2/2 € A m and 
T < t\ < <2 < T' satisfying the following Chapman-Kolmogorov axioms: 

(CK1) U m {y 1 ,t 1 ;y 2 ,t 2 )>0 Vj/i,t/ 2 eA, Vtj < t 2 € M, 

(CK2) U m (y 1 ,t 1 ;y 2 ,t 1 ) = S yiy2 Vyi,y 2 GA, Vt x G M, 

(CK3) ^2 U m (yi,ti;y2,t2)U m (y2 ) t 2 ;y3,t 3 ) = U m (yxM)V3M) 

Vyi,y 3 e A, Vii < t 2 < t 3 S E. 

Given a Markov propagator U m (yi,ti;y 2 ,t 2 ), one can define a sequence of n— point function 
having all the necessary properties by setting 

(5.1) $m(yi,t 1 ;...;y n ,t n \y )= ]J U^yj-i^j-^yj,^). 

j=l..n 

where t = T. 

Definition 19. (Markov Process.) A filtered probability space generated by a Markov propa- 
gator is called Markovian and the corresponding stochastic process is called Markov process. 

Definition 20. (Markov Generator.) If the matrix elements U m (yi,ti;y 2 ,t 2 ) are differen- 
tiable functions of the time parameter t 2 in a right neighborhood of t\, then one defines the 
Markov generator at time t\ as the following right derivative: 

(5.2) £ ro (yi,y 2 ;£i) = lip -^-U m (yi,tr, y-2,h) 

tzltt dt 2 

Proposition 1. // C m (y\, y 2 ; ti) is a Markov generator, then for all pairs 2/1,2/2 € A m the 
following two properties hold: 

(5.3) (MG1) C m ( yi ,y 2 ;t)>0 if y x + y 2 , 

(5.4) (MG2) C m ( yi , yi ;t) =- £ C m (y u y 2 ;t). 

Viceversa, if C m {yi,y 2 ;t) is a differentiable one-parameter family of matrices satisfying condi- 
tions {A) and (B) above, then the differential equation (5.2) admits one and only one solution 
satisfying the initial condition U m {yi,t;y 2 ,t) = S yiy2 . 



The propagator U m (yi,ti;y 2 ,t 2 ) defined by the differential equation in (5.2 1 can be repre- 
sented by means of a so-called path- exponential defined as follows. Let N > be an integer and 
let us consider the product 



(5.5) U^(y 1 ,t 1 ;y 2 ,t 2 )= ( l + ^£ ro (t x ) I ( 1 + 8t N C m (t N ) 
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where Stiy — ^jy 1 ■ If AT is so large that 

(5.6) (6t) N < min C m (y, y; 

i/eA,te[*i,t 2 ] 



then the operator product U^(yi, t\\ y 2 , t 2 ) in equation (5.51 is a probability kernel. Passing to 
the limit N — > oo, we find 

(5-7) Jim U%(y 1 ,t 1 ;y 2 ,t 2 ) = U m (y 1 ,t 1 ;y 2 ,t 2 ). 

By expanding the product in ( |5.5[ ) and passing to the limit, we arrive at the following: 

Theorem 3. (Dyson expansion.) The probability kernel can be represented as the following 
convergent series: 
(5.8) ' 



U m (yi,h;y 2 ,t 2 ) = (l + ^2 / ds i / ds 2 ... ds n C m (s{)C ri 



,(s 2 )...X m (s„) (yi,y 2 )- 



By differentiating with respect to the two time coordinates in the Dyson expansion, we also 
find the following two equations: 

Theorem 4. (Forward and backward equations.) The probability kernel satisfies the back- 
ward equation 

(5.9) J-t/m(*i;*2)+An(*l)E/m(*i;*2) = 

as well as the forward equation 

d 

(5.10) ^U m (h;t 2 ) = U m (h;t 2 )C m (t 2 ). 

ot 2 

A handy notation for this expansion is given by the following: 



Definition 21. (Path-ordered Exponential.) The equation (5.8) is written as a path ordered 
exponential 

(5.11) U m {y x ,ti\y 2 ,t 2 ) = Pexp( / £ m (s)ds)(yl,y2) 



where the operator P formally acts as follows: 



(5.12) P / C m (s)ds) =n\ dsi / ds 2 ... / ds„C m (si)C 

m (s 2 ) £ m (-S n ) . 

Wti / Jti JS! Js n -i 

Using the Dyson expansion for the path-ordered exponential, one finds a path-integral repre- 
sentation of the probability kernel. Let us set the following definition: 

Definition 22. (Symbolic Path.) A symbolic path 7 = {70,71,72, ■•■•} is an infinite sequence 
of sites in h m 1j such that 7^ ^ 7y_i for all j — 1, .... Let T m be the set of all symbolic paths in 
h m 1. 

Theorem 5. (Path-Integral Representation.) The propagator admits the following repre- 
sentation: 
(5.13) 

U m (x,T;y,T') ^2 ds i ds 2 ... ds q 

q=l -Y&r m :y =x,y q =y Jt J Sl J s i~ 1 

(5.14) 6XP (/ £m (T 'T°' ,,|1 )* 1) j II' C ' r, ^- 1, ^ ;s ^ exp (/ £m(lj,lj;vj)dv. 

where t q+ i = T . 
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Notice that the total mass of the sector of path space H^yo, T; T") consisting of lattice paths 
originating from y <E A at time T and attaining at most q different values by time T' is given 

by 



KWL(yo,T;T')) = V [ ds x [ ds 2 ... [ ds q exp([ Cmilo, 7o; v )dv ) 
7 er m:70 =x- /T J 'i Js «-i ^ Jt ' 

IJ £ ™(7j-i,7j; *j) exp / Cmi-yj^^v^dvj . 

j=l \ Js i ' 

where s m+1 = T". If A m = max :Eej 4 7ii |£ m (x,x)| and c = sup te [ TT /]||£ TO (f)|| is the operator norm 
of the matrix £ m (t), then we have that 

(5.15) n[Hl(y ,T;T>)] < ^ ~ T ) q e -^(T>-T) 

This expression reaches a maximum as a function of q at q w c(T" —T) and then declines at super- 
exponential rate. Hence, by discretizing the space coordinate, we ensure that the probability 
mass of the set of paths with q changes decreases faster than any exponential in the limit as 
q — > oo. This convergence is however not uniform as m — > oo and h m — > 0. 

Definition 23. (Inverse lattice.) Le£ us consider the following lattice: 



(5.16) 5 m =j z h — = 0, ..2 m - 1| 

a/so ca//ed Brillouin zone or inverse lattice with respect to A m . 

Definition 24. (Pseudo-differential Symbols.) The symbol of a Markov generator is defined 
as follows: 

(5.17) £ m (x,p;t) = £ m (x,y)e ip ^ 

yeA m 

where p € B m . 

Two particularly important special examples of Markov processes are given by monotonic 
processes and diffusions. 

Definition 25. (Monotonic process.) A Markov process of generator £ m (x,y) on the sim- 
plicial sequence A m is said monotonic non- decreasing if 

(5.18) £ m (x,y)=0 whenever y < x. 
and monotonic non-increasing if 

(5.19) £ m (x,y)=0 whenever y > x. 

Definition 26. (Diffusion Process.) A diffusion process is a Markov process with generator 
of the form 

(5.20) £ m (x, y; t) = » m (x; t)V hm (x, y) + A hfn (x, y) 
where 

(5.21) V h (x,y)= and Afe(a . y) = W + &v£-h ~ ™ y , x 

and fi(x;t) — (fi m (x;t)) and cr(x;t) = (a m (x;t)) are two simplicial functions which, for simplic- 
ity, we assume smooth in both arguments. 

Proposition 2. (i) The symbol of a diffusion process is given by 

~ s'mph cr m (x;t) 2 cosph - 1 
(5-22) £ m (x,p;t) = (J, m {x;t)— — + ~ 2 . 
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(ii) The path integral representation for a diffusion process has the form 

oo 

(5.23) U m (x,T;y,T') =^~ q E W ro ( 7 ,g;T,T') 

9=1 7 € r m : 70 = z,7 g = ?y 
|7i-7 J --i|=l Vj>l 

w/iere 

W m { 1 ,q;T,T')= ( d Sl [ ds 2 ... [ ds q 

JT Jsi Js q -1 

(5.24) X c <"<™^*»> JJ / £ -^^^)^2£ ro (7 j ,7 j+1 )) 
and so = 7 1 - 

6. Martingales and Monotonic Processes 

In this section we introduce the notion of piecewise smooth Markov process which covers 
a large family of models useful for applications. In this context, we define martingales and 
monotonic Markov processes. 

Definition 27. (Piecewise Smooth Markov Processes.) Consider the time interval [T,T'j, 
the simplicial sequence A m , m > m and a finite number of time points T = t < t\ < .. < t n = 
T'. A piecewise smooth Markov process is given by a family of Markov generators C % m (y\,y2\t) 
defined on each half open time interval where i = 0,1, ...n— 1. In correspondence to 

each time point ti,i = 1,2, ...n, one also defines a mapping operator U^ n (x,y) such that 

(6.1) (mai) w m (yi,y2)>o 

(6.2) (MA2) J2 u L(yi,y2) = i v yi g A. 

J/2 

The Markov propagator for any pair of time points ti < s < s' < t i+1 is defined as follows: 
(6-3) U m (y 1 ,s;y 2 ,s') = Pexp ^ £ l m (v)dvj {yi,y 2 )- 

Moreover, if s' = t i+ i, then 

(6.4) U m (y 1 ,s;y 2 ,t i+1 ) = ^ Pexp ^ C m (v)dv^j (2/1,2/3)^(2/3,2/2)- 



ya6A 

More general Markov propagators are obtained by taking products of the ones above. 

Definition 28. (Attainable Sets.) Let U m be a piecewise smooth Markov process, m > m , 
y G A m and t G [T, T'\. The attainable set T> m {U, t, y) C A m is defined as follows: ift G (U, ti+i) 
for some i = 0, ..n — 1, then T> m (U, t, y) is the set of y G A such that Li{y, y; t) > 0. // instead 
t = ti for some 2 = 1, ..n — 1, then V m (U,ti,y) is defined as the set of the y G A such that 
Ui{y,y;t) > 0. 

Definition 29. (Equivalent Markov Processes.) Two piecewise smooth Markov propagators 
U and U' are called equivalent if their attainable sets T> m (U,t,y) and T> m (U' ,t,y) are equal for 
all t G [T,T'\ and all y G A. IfV m (U,t,y) is a subset ofV m (U',t,y) for all t G [T,T'\ and all 
y G A, then one says that U m is absolutely continuous with respect to U' m . 

Definition 30. (Measure Changes.) Let U m ,m > mo be a family of piecewise smooth Markov 
propagators. A measure change is characterized by a family of positive, non-zero functions 
Gm(y') — indexed by y G A m and t G [T,T'] t which is strictly positive for all y' G T> m (U,y,t) 
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and zero otherwise. A measure change junction defines a transformation of a Markov generator 
into an equivalent one according to the following formula: 

(6.5) C'(y,y';t) = -i_£( y , y ' ; ^0^(2/') - -*(£(t)G*)( y ) 5 yy ,. 

Notice that the specification of the function G v j l (y') at the point y' — y is immaterial in the 
sense that it does not affect the measure change transformation. 

Definition 31. (Time Changes.) A measure change is called time change if there is a function 
4>{y,t) such that G^(y') = 4>(y,t) for all y G A m , y' G T> m (U,y,t) and t G [T,T']. The time 
change is called state independent if 4>(y, t) = 4>{t) is a function of the time coordinate only. 

Theorem 6. (Deterministic Time Changes.) If (j>(t) defines a state independent time change 
corresponding to the measure change function G^*(y') = <fr(t), then 

(6.6) U' m {y, t; y', t') = U m (y, X(t);y', X(t')) 
where 



(6.7) X(t) = J cf>(s)ds. 

The following is a particularly interesting special case of measure change: 

Definition 32. (Numeraire Changes.) Consider a smooth (as opposed to just piece-wise 
smooth) Markov process of generator C m (y; t). Let g m (y'', t) be a function satisfying the equation 

(6.8) d J^l + {Cm g m ){ m t)=Q. 

The measure change given by the function G^(y') — g m (y';t) is called numeraire change and 
the Markov generator transforms as follows: 

ran\ r> ( ' *\ 9m(l/\t) , , . 1 dg m (y;t) 
(6-9) £ m (y,y;t) = -—-£ m (y,y;t) + — — — 5 vy ,. 

g m {y;t) g m {y\t) dt 

Notice that if g m {i) denotes the multiplication operator of kernel g m {y\ t)8 yy i, then equation 
(13.41 can be written more compactly as follows: 

(6.10) C' m {t) = -J_£ m (% ro(t ) + 1 d g " (t) . 

g m {tj g m \tj at 



Theorem 7. (Numeraire Changes.) If g m satisfies equation (13.3) and defines a numeraire 
change, then 

(6.H) U' m {y,My\t') = 9mi f^ U m (y,t;y'^). 

9m{y,t) 

Let Ft — F(y.,t) be an adapted process in the time interval [T,T']. Let us consider a fixed 
time t G [T, T"). Recall that, since F t is non-anticipatory we have that F(y.,t) — F(y',t) for 
all pairs of paths such that y s — y' s for all s < t. If y £ T> m (U' ,t,y), let y. = Extt(y. ,yt) be 
the constant extension path such that y s = y s for all s < t and y s — yt for all s > t. With 
probability one, we have that 

(6.12) Urn F{y..t + St) = F(Ext t {y.,y),t) 

ot[0 

for some y G T> m (U, t, y). 

Definition 33. (Monotonic Processes.) Let U be a piecewise smooth Markov propagator and 
let F t be an adapted process given by the non- anticipatory path functional F(y.,t). F t is said to 
be increasing at time t if 

(i) For all y G V m {U,t,y) we have that 

(6.13) F(Vxt t (y.,y),t)-F(y,t)>0. 
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(6.14) 



(ii) We have that 

dF(Ext a (y.,y t - ),s) 



> 



ds 

where y t -o = limstioyt-St- 
(iii) For all y G A m and all t G [T,T'], either the inequality in (i) holds in a strict sense for 
at least one y G T> m (U,t,y) or the inequality in (ii) holds in a strict sense. 
In case the property (iii) fails but the other two still hold, the process is called non decreasing. 

Theorem 8. (Monotonic Processes.) Let U m and U' m be two piecewise smooth Markov prop- 
agators and let F{y. , t) be a non- anticipatory path functional. If U m and U' m are equivalent and 
if F(y. , t) regarded as an adapted process under the path measure generated by U m is increas- 
ing (non-decreasing), then also F(y.,t) regarded as an adapted process under U' m is increasing 
(non- decreasing) . 

Proof. This theorem descends from the fact that the definition of monotonicity depends only on 
the attainable sets. □ 

Definition 34. (Martingale Processes.) A process F t is a martingale if for all times t we 
have that 

(6.15) J^(Efart t (tf.,ite_ ),t)+ £(yt,y,t){F(Eid t (y.,y),t) - F(Ext t (y.,y t _ ),t)) =0 

VBV(U,t,y t ) 

The process F t is called an equivalent martingale if there is a second piecewise smooth Markov 
propagator U' for which the non- anticipatory path functional F(y. ,t) is a martingale process. 

Theorem 9. (Equivalent Martingales.) 

(i) If F t is an increasing adapted process then it is not an equivalent martingale. Otherwise 
stated, if F t is an equivalent martingale then it is not increasing. 

(ii) If F t is a non-decreasing adapted process and it is also an equivalent martingale, then it 
is a constant process with F(y.,t) — const for all t G [T,T']. 

Proof. This theorem descends from the fact that monotonicity properties are preserved by mea- 
sure changes. □ 

Martingales are particularly useful as they can be constructed by taking expectations. 
Let <&(?/.) be a continuous path- functional. From the modeling viewpoint, such a functional 
can represent future cash flow streams. For instance, one choice could be 

(6.16) *(y.) = *o(Wf) 

where $o is a continuous univariate function and t' G [T, T'\ is fixed. In a more general example, 
one may consider a path functional of the form 

(6.17) $(y.) = / ds 1 ... ds n a(t;si,....s n )F 1 (y Sl ,si)...F n (y Sn ,s n ) 

with a(t; si, ....s n ) = for t G [T, i']. The path conditioned expectation of ^(y.) is the non- 
anticipatory path functional F t — F{y., t) such that 



(6.18) F(y.,t)= / *{z.)n[dz.\ 

JA(y.,t) 

where the integral is restricted to the set A(y.,t) = {z.\z s = y s Vs < t}. The intersection of 
the set A(y. ,t) with each of the spaces r H.% m is a compact, finite dimensional submanifold with 
boundaries. To denote path conditioning, we also use the following notation: 

(6.19) F t = F{y.,t) = E[*{y.)\y» s<t] = E t My.)}. 



Proposition 3. The process F t in {6.19) is a martingale for t <t\. 



Proof. Equation (6.15) descends from the backward equation in (5.9 ) . □ 
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7. The Fundamental Theorem of Finance 

In this section, we derive the Fundamental Theorem of Finance in a context general enough 
to encompass most cases of practical relevance. 

Definition 35. (Financial Model.) Let A m ,m > mo, be a simplicial sequence. A financial 
model is given by a family A^(y.,t), k = 1, ...n of adapted processes modeling asset prices and 
a non-decreasing adapted process B m {y.,t) modeling the money-market account. For notational 
convenience, we set A^(y.,t) — B m (y.,t). Let us introduce also the discounted asset price process 
defined as follows: 

A k m (y.,t) 



(7-1) At{y.,t) 



B m (y-,t) ' 



Definition 36. (Trading Strategies.) Given a financial model with n assets, a trading strategy 
is given by a family of adapted processes (!^(y.,t), k — 0, ....n. The value process of a strategy is 
the adapted process 

n 

(7-2) U m (y.,t) = J2C(y;t)A k m (y.,t) 

k=0 

The discounted value process instead is given by 

n 

(7.3) u m ( y .,t) = J2C(y-M k m (y;t) 

k=0 

Definition 37. (Self-Financing Condition.) An adapted trading strategy is called self- 
financing if the following two conditions hold: 

(7.4) (SFl) J2^(Ext t (y., yt ^),t)A k n (E X t t (y.,y t ^),t) = 

k=0 

and if, for all y € T> m (U, t,y t ), we have that 

n 

(7.5) (SF2) 2(C*(E X t t (y.,y) ) t)-C*(Ext t (». ) y t _o),t))^(E x t t (y.,»),t) = 0. 

k=0 

Proposition 4. IfC!^(y.,t) is a self-financing trading strategy, then the corresponding discounted 
value process satisfies 

(ij a\ dn m (Ext t (y.,2/t-o),*) VVfc /pv + „ ^ t\ ( Ext t fa- , yt-o),t) 
(7-D) Q~ t = 2^Cm(Ext t (y.,?/t_o),t) ^ 

fc=0 

and for all y € T> m (U, t, y t ), we have that 
(7.7) n m (Ext t (j/.,j/),t) -n m (Ext t (y.,y t _ ),t) 

n 

(7-8) = Cm (Ext t (y. ,y),t) (A k n (Ext t (y. ,y),t) - A k m (Ext t (y. ,y t . ),t)). 

fc=0 

Definition 38. (Arbitrage Strategies.) A self- financing strategy is called arbitrage at time 
t if the corresponding discounted value process Il m (y.,t) is increasing at time t. 

Theorem 10. (Fundamental Theorem of Finance.) If there is an equivalent measure with 
respect to which all discounted base asset price processes are martingales, then 

(i) The discounted value process of any self- financing trading strategy under the same equiv- 
alent measure is a martingale. 

(ii) There is no arbitrage. 

Conversely, if there is no arbitrage than there exists a measure change under which all discounted 
asset price processes become martingales. 
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Proof. The first part of the theorem is a simple consequence of the definitions put forward. The 
converse instead requires a proof. 

Let V m (U, t, yt) be the vector space spanned by functions of the form v(y) where y £ T> m (U, t, yt)- 
Let K be the cone in V m (U,t,yt) made up by all the vectors with non-negative components. 
Also, let us introduce the following vector £^ £ V m (U,t,yt): 
(7-9) 

C(V) = Sy, yt _ 9J ^ {Extt ^ ,yt - o),t) + (1 - Sy, m _ )(A k m (Ext t (y.,y),t) - A k m (Ext t (y.,y t . ),t)). 

Notice that, in case k = and the asset is the money market account, then ^(y) = for all 
y£V m (U,t,y t ). 

Suppose there is no arbitrage and fix a time t. Assuming there does not exist a trading 
strategy with a strictly increasing value process, for all vectors (v fc )fc=i....n £ V m (U,t,y t ) there 
are two elements y+,y_ £ T> m (U, t, y t ) such that 

n n 

(7.10) ]>> fc ^(y + )>0 while 5> fe £*(y_) < 0. 

k=l k=l 

If in equation (7.101 we set v k = Ski, we conclude that the vector {£,} n {y))y £V njtyA ^ as 
both positive and negative components. Hence, the the hyperplane II 1 C V m (U, t, y t ) orthogonal 
to the vector £^(y) £ V m (U,t,y t ) intersects K. 

Let Pi be the orthogonal projection operator onto the hyperplane II 1 and let us consider the 
vector 

(7-11) (-PlO(y) = £m\V) - v=^ -PL r=\ tm(y)- 

2-/zeV m (U,t,y t ) im\ Z ) l im\ Z ) 

for i = 2, ...n. Due to absence of arbitrage, there are two elements £ D m (U,t,y t ) such 

that 

(7.12) (P 1 e m )(y+)>0 while {P 1 e m ){y-) <Q- 

Hence, the vector Pi^ is transversal to the octant K of positive vectors. As a consequence, the 
hyperplane U2 C HiV is orthogonal to both vectors and also intersects K. 

The argument above can be iterated n times, leading to the conclusion that there exists a 
strictly positive function gl£(y) > 0, y £ V m (U, t, y t ), such that 

(7.13) ]T e m (y)gt(v) - 

yev m (u.t, yt ) 

for all k = 0, ...n. In particular, this implies that there exists a measure change function G^(y) 
such that 

(7.14) Y, ( Ah m(Vxt t (y. ,y),t)~ A^(Ext t (y., y t -o),t))jC(y, y; t)G£{g) = 0. 

yehl. d 

□ 



8. Weak Convergence of Markov Generators 

Consider a one dimensional Markov process defined on the simplicial sequence A m , m > mo 
and the time interval [T, T'] . Many different specifications of Markov generators on A m may 
correspond to the same limit. In this section we identify a canonical sequence of generators 
under a few regularity hypotheses which imply the existence of a weak limit in distribution sense 
for the generator. This is a necessary first step to single out the general form of an admissible 
Markov generator. In the following sections, we then investigate convergence under the much 
finer criteria of pointwise convergence for probability kernels and their derivatives. 

First consider the case when the limiting domain Aao is bounded. Without restricting gener- 
ality, let us suppose that Aoa — [—L,L]. 



OPERATOR METHODS, ABELIAN PROCESSES AND DYNAMIC CONDITIONING 



17 



Let L m be a sequence of Markov generators. The first assumption we make is that the first 
two moments are finite, or more specifically 

Hypothesis MG1. The sequences 

(8.1) n m (x,t)= ^2 £m{x,y;t)(y - x), 

yeA m 

are uniformly bounded in absolute value as m 
all t € [T,T'], i.e. the following limits exist: 

(8.2) [i(x, t) = lim fi m (x,t), 

m— *oo 

Notice that, due to the dominated convergence theorem, the functions /j, m (x,t) and a m (x,t) 
regarded as piecewise constant functions on A x converge weakly to the corresponding limits. 

Hypothesis MG2. The following limits exist and are finite for all m > m and all pairs 
x, y e A m such that \x — y\> 2h m : 

(8.3) X m (x,y;t)= lim V £ m ,(x,z;t). 

m' — >oo * — ' 

zEA m / :y — h rn <z<y-\-h rn 



<r m (x, t)= / X] ^{x, y'' ~ x ) 2 

V yeA m 

oo and converge to limits for all x e A m and 
a(x,t) = lim a m (x,t). 



The family of functions A m (x, y) can be represented in terms of the so called Levy measures 
specified as follows: 

Theorem 11. (Levy measures.) For all x e A m there exists a measure v xt {d(;) in [— L, L] 
(called Levy measure) with the following two properties: 

(i) 

(8.4) J^ xt (doe<™ 

(ii) For all y G A m with \x — y\> 2h m the following representation is valid: 

ry+h m /2 

(8.5) \ m (x,y)=Km v x t{dQ- 

£ J-° Jy-h m /2+E 

Definition 39. (Finite activity jumps.) A Markov process is said to have finite activity 
jumps if for all m > mo and all x e A m we have that 

(8.6) J i/ xt (d^)<oo 

Notice that given MG1 and MG2 we have 

(8.7) limsup^2 \ m (x,y;t)(y - x) 2 < C m (x, y; t)(y - x) 2 < oo. 

m — >oo * 

y y 
More generally, if e V(R) is a test function such that 0(0) = and 4>'(0) = 0, then we have 
that 

(8.8) limsup^ A m (x, y; i)<j>(y — x) < oo. 

m — >oo 

y 

Although the above sequence is bounded, the limit may not exist in general. We thus need to 
stipulate this as a separate assumption: 

Hypothesis MG3. For all test functions <fi € T>(R) such that </>(0) = (f>'(0) = and all x e A m , 

the sequence ^2 X m >(x, y; t)<j)(y — x) admits a limit as m! — > oo. 
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Let us introduce the following notations: 

(8.9) l&(x-,t)=Y i X fn (x,ir,t)to-x), a° m (x;t) = ^\ m (x,y-t){y - xf . 

Notice that the sequence a^ n {x,t) is non-negative and is uniformly bounded as a function of 
to. More precisely 

(8.10) <a° m (x,t) <a m (x,t). 

On the other hand, we cannot conclude in general that ^(i, t) is necessarily uniformly bounded. 
In fact, the difference ii m {x,t) — /^(x,i) could diverge as m — > oo while being compensated by 
terms concentrated at y — x ± 1 which in turn also diverge in the limit as m — > oo while keeping 
the total drift fi m (x,t) uniformly bounded. An exception to this general situation is found in 
case the following condition is satisfied: 

Hypothesis MG4. For all m > m Q and all x E A m we have that 

(8.11) J \x\v xt {d£) <oo 



Two particularly important situations in which this condition MGA holds are given in the 
following: 

Theorem 12. (MG4.) Assuming MG1, MG2, MG3, if the Markov process is either monotonic 
or has finite activity jumps, then MGA holds. 

Under the three assumptions MG1, MG2, MG3, the sequence C m of Markov generators can 
be mapped into an equivalent canonical sequence of Markov generators C . More precisely, we 
set 

(8.12) £C{x,y;t) = \(x,y;t) 
in case \x — y\> 2h rn . Furthermore, we set 

(8.13) Cg(x, x ± h m ; t) = jl m (x; t)V hm (x, x±l)+ A hm (x, x±l) 
where 

(8.14) Ji m {x;t) = n m {x) - pp m {x), a m (x;t) 2 = a m (x;t) 2 - <j° m (x;t) 2 . 
and the operators and A^ are defined as in equation ( |5.21 1. Finally, 

(8.15) Cg{x,x;t) = Yl 

Theorem 13. (Canonical Representations of Markov Generators.) For all smooth func- 
tions of compact support <f> S V^A^) and all x £ A m we have that 

(8.16) lim V CF m , {x,y)<P{y)= lim V C m ,{x,y)<l>{y). 

m — >oo L — ' m' — >oo L — ' 

A compact representation for a canonical generator summarizing these definition is obtained 
by considering the symbol as specified in Definition ( 24 1 . 

Theorem 14. (Levy-Khintchine representation.) Under the hypothesis MG1, MG2, MG3 
above and in case A^ — [-L, L] is bounded, the symbol of a generator in canonical form can be 
expressed as follows: 
(8.17) 

Cg(x, P ; t) = i^x; 1) S ^+a(x; tf C ^^ + £ (e*»-) - 1 - i^(y - x)) X m (x, y; t). 
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The following limit converges in the weak topology: 

(8.18) lim tZ{x,p-t) = ip{x-t)p-^^p 2 + [ (e<*-l-tpO"xt(dO> 

m— >oo A J — L 

where v x t{d£) is the Levy measure. Moreover, if also MGA is satisfied, then the symbol of a 
generator in canonical form can be expressed as follows: 

(8.19) cZ^P;t)= i Kx;t)^ + a(x;tr C0SP ^- 1 + £ (e^y-^-l)\ m {x,y;t). 
The following limit converges in the weak topology: 

(8.20) lim £f n (x,p;t) = ift(x; t)p - ^^-p* + C ( e ** - l)p xt {<®- 

m— >oo 2 J — L 

9. Fast Exponentiation and Spectral Methods 

Numerical analysis of pricing models in the operator formalism depend on the ability to com- 
pute the propagator for a given generator C(t). Time homogenous Markov generators C m (y\, y 2 ) 
represent a privileged special case of particular importance. In this case, the associated propa- 
gator solving the differential equation in ( |5.2[ ) is given by the matrix exponential 

(9.1) U m (y 1 ,s;y 2 ,s') = exp((s' - s)£ m )(y 1 ,y 2 )- 

Matrix exponentiation can be defined in several equivalent ways, such as by Taylor expansion 

00 P 

(9.2) exp(i£ m ) = j^L 
or by means of Neper's formula 

(9.3) exp(i£ m ) = lim (l + ^-C r , 

N^oo \ JV 



N 



We find that the most efficient and robust method for exponentiating Markov generators is 
the so called fast exponentiation algorithm. Let us fix a Markov generator h{x,y) and a time 
horizon t. Let 5t > be the largest time interval for which both of the following properties are 
satisfied: 

(FE1) min(l + 8tL(y, y)) > 1 12 

(FE1) log 2 ^=nel 

To compute e tc (x,y), we first define the elementary propagator 
(9.4) ust(x,y) = S xy + 5tL(x,y) 

and then evaluate in sequence u 2 st = ug t ■ Ust, u^t = u 2 st ■ u 2St , ... u 2 n St = u 2 ^-i St ■ u 2 n-i st . 

As we show in the next few sections, this algorithm approximates probability kernels with 
errors with respect to the continuous time kernel density which, in fairly general cases, are of 
order 0(ti^ n \\og h m \ 3 ). This is the same order by which the continuous time kernel density differs 
from its continuous limit also according to estimates below. In the case of Brownian motion, a 
sharp convergence estimates of ordered 0(/i^J can also be proven. 

Matrix multiplication is accomplished numerically by invoking the routine dgemm in Level- 3 
BLAS. Very efficient, processor specific version of dgemm are now available along with implemen- 
tations on massively parallel GPU chipsets. It turns out that the standard measure of algorithmic 
complexity as the number of floating point operations required to accomplish a certain task, is 
not simply proportional and scales non-linearly with respect to execution time. Using blocking 
and cache optimizations and distributing the load across many cores, execution time for medium 
to large matrices appears to scales much better than the naive n 3 scaling one would obtain by 
triple looping (Goto and van de Geijn to appear) . 
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A more general method that allows to compute not only exponentials but also other func- 
tions of a Markov generator is full diagonalization. Unfortunately, unless the Markov generator 
is symmetrizable, diagonalization algorithms can possibly run into serious instabilities within 
double precision arithmetics because of the phenomenon of pseudo-spectrum (Trefethen and 
Embree 2006). This makes it impossible to diagonalizc exactly within the limits of double- 
precision arithmetics and forces one to resort to expedients such as for instance spectral trun- 
cations. Since the fast exponentiation method has much better scaling properties than full 
diagonalization and is entirely stable, we recommend that it be used in all situations where a 
matrix exponentiation is required. However, it often arises the necessity to define and com- 
pute also other functions of a given time independent Markov generator and for this purpose 
diagonalization can be usefully employed, especially when the target matrix is symmetrizable. 

Not all Markov generators are diagonalizable but most are and it is safe to take diagonaliz- 
ability for granted in numerical applications. To make this statement more precise, one sets out 
the following definitions: 

Definition 40. (Generic Properties.) A dense G$ of a topological space is a countable 
intersection of dense open subsets. If a property is valid on a dense Gs one says that it is valid 
generically, or that it is generic. Instead, if the same topological set is endowed of a measure 
and a property is valid on a full measure set of parameters, one says that it is valid almost surely 
with respect to that particular measure. 

We have that 

Proposition 5. Markov generators are diagonalizable both generically and almost surely. 
Let u n (x) and A„, n = 1, ..N, be eigenfunctions and eigenvalues of the Markov generator C, 

i.e. 

(9.5) Cu n = X n u n . 

Let U be the matrix whose columns are given by the vectors u n (x) and let A be the diagonal 
matrix with the eigenvalues A„ on the diagonal. Hence 

(9.6) CU = UA. 

We have that C = UAU -1 and e tc — Ue tA U~ 1 . This equation expressed in components reads 
as follows: 

N 

(9.7) e tc (x, x') = e Xnt u n (x)v n (y) 

71=1 

where v n (y) is the n— th row vector of the inverse matrix V = U^ 1 . 

One may extend the above definition to other functions of a Markov generator. If tp(X) is a 
function, one may define "4>(£) as the operator whose matrix is given by 

N 

(9.8) ^(C)(x,x') = ^(X 

It 

n=l 

An important example of functional calculus is found to express the Markov generator of pro- 
cesses obtained by stochastic time change, whereby the time-change process has independent, 
uniformly distributed increments. Processes in this class are called Bochner subordinators. Be- 
cause of time and space homogeneity, Bochner subordinators can be constructed starting from 
a process on simplicial sequence /i m Z and are characterized by a Markov generator of the form 

£ m {x,y) = £h m {y - x). 

Theorem 15. (Bochner Subordinators.) If the limit lim/^o^/i = I exists in weak sense on 
V(R) then the limit kernel has the following form 



(9.9) J £{x)<j>{x)dx = n</)'(0) + J 



(<f>(x) - mHdx) 
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where v{dx) is a positive measure supported on M. + and is such that 

POO 

(9.10) / xv(dx) < oo. 



The characteristic function of the process defined as the Fourier transform of the Markov 
generator is thus 

/•oo 

(9.11) e(^limVe rtl 4W = #+ (e lkx - l)v{dx). 

h i-° x JO 

The Laplace transform instead is given by 

/•OO 

(9.12) (f)(\) = -e(i\) = -lim^2e- Xx l h (x)= n\+ {I - er Xx )v{dx). 

o x Jo 

A function of this form is called Bernstein function. 

Bernstein functions may be used to express Laplace transform of the transition probability 
kernel of time-homogeneous monotonic processes as follows: 

poo 

(9.13) / e tc (0,x)e' Xx dx^e^ x \ 



In turn, we may also express the transition probability kernel in terms of the characteristic 
function, i.e. 

f°° rib 

(9.14) e tc (0,x) = / e lte « +lfea: r 



— CO 



2tt 



Notice that since </>(— ik) — —e(k), the last formula may also be interpreted as the Fourier-Mcllin 
inversion of the previous one. 

Example 1. (Poisson Process) The Poisson process corresponds to 

poo 

(9.15) 4>p(X; c) = c(l - e~ A ) = c / (1 - e - xt )6(t - \)dt. 

Jo 

Example 2. (Stable Process) The stable subordinator with index a € (0, 1) is given by 



(9.16) 0s (A;a) = A a = — -/ (1 - e- At )r 1_a *- 

r U -") Jo 

Example 3. (Gamma Process) The Gamma subordinator with variance rate v > is given 
by 

1 1 c°° 

(9.17) <WA;z/) = -log(l + i/A) = - / (1 - e-^t^e^^dt. 



To add jumps to a diffusion process one can use the method of independent stochastic subor- 
dination. Namely let F t be the diffusion process with drift function fJ,(F) and volatility function 
<y{F) and let L (?/i, y2', cr) be the corresponding Markov generator. Consider the generator 

(9.18) Lv G (/i, cr, v) = -</v G (-L d (/i, a); v). 

The propagator for L J (/i, cr, v) satisfies the following equation: 
(9.19) 

JV 

exp(tL VG ( f i,o-,nu))(y 1 ,y 2 ) = E [exp (TfL d ( f i,cr))(y 1 ,y 2 )] = £ e^^^^u^v^y) 

n=l 

where T" are the paths of the monotonic process in Example 3, the u n {x) arc the cigcnfunctions 
of the diffusion generator L d (/x, cr), A„ are the corresponding eigenvalues and the functions v n {x) 
are defined as in (9.7l. Hence, the generator Ly G (/i, a, v) identifies the time-changed process 
with paths Ft». 
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To model asymmetric jumps one can follow several strategies. A simple one is to specify the 
two different variance rates v + and v_ for the up and down jumps and compute separately two 
Markov generators 

(9.20) L VG ( M , o,v±) = -<f> VG (-L d (n, a);v±). 

The new generator for our process with asymmetric jumps is obtained by combining the two 
generators above 

{ LyG (j/i, 3/25 if F(y 1 ) < F(y 2 ) 
L V G(2/i,ZV2;/i,CT,i/ + ) if F(y x ) > F(y 2 ) 

The construction can also be localized. If <j) x (X) is a family of Bernstein functions indexed by 
the state variable x G A m , then one can consider the Markov generator of matrix 

(9.22) l(x, y) = —cj) x (—h d (fj,, a))(x, y). 

This construction allows one to model state dependent jumps. 

10. Construction of Brownian Motion 

This section is based on work in collaboration with Alex Mijatovic, see (Albanese and Mijatovic 
2006). 

Let A = (A m ),m = to , to +1, ... be a simplicial sequence converging to an interval lim m ^oo A m 
[-L, L] C R, where < L < 00. Let \i and a be two constants. 
The generator of a Brownian motion on A m has the form 

(10-1) An - MV ftm + \<J 2 A hm 

for all interior points x G Int(A m ). We assume that m is large enough so that 

(10-2) £ > lf " 



2d?,, 2A m 
for all m > too- 

If a; is a boundary point, i.e. x € d(A m ), then several definitions of the operator C m are possi- 
ble depending on the choice of boundary conditions. Absorbing boundary conditions correspond 
to the choice 

(10.3) £ m (x,y) = VyeA m . 
Reflecting boundary conditions correspond to 

(10.4) C m (x,x) = Li\7 hm (x,x) + ^o- 2 A hm (x,x) 

(10.5) C m (x, y) = -C m (x, x) 

where y G Int(A m ) is the closest point to x in the interior of A m , while C m {x,y) = for all 
other points y. Periodic boundary conditions are implemented by setting 

(10-6) £ m (x,y) = ii\7 hm (x,y) + ^cr 2 A hm (x,y) 

if y = x or y G Int{A m ) is the closest point to x in the interior of A m , and also 

(10.7) C m {x, x 1 ) = nV hm (x, y) + \o 2 A hm {x, y) 

where x' G d(A m ) is the boundary point at the opposite extreme of the simplex A m . Finally, 
mixed boundary conditions can also be defined by taking a convex linear combination of the 
generators satisfying one of the three boundary conditions above. 
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Theorem 16. (Convergence Estimates for Brownian Motion.) Consider a Brownian 
motion defined as the process with generator in (10.1) where [i and a are constants satisfying the 
bounds in (10.2) for all m > ttiq and we assume periodic boundary conditions. Let us consider 
the kernels 



U m (x,0;y,T) = e m (x,y) 



(10.8) 

and 

(10.9) U%(x,0;y,T) = (l + StcJ)^ (x,y). 

1 h 2 

If 5t is such that < 5t < | 2a 2 -h\p,\ an< ^ m ° * s lo,rge enough, then there is a constant c > such 
that for all m! > m > mo we have the following inequalities: 

(i) 

(10.10) 

(ii) 
(10.11) 
(iii) 
(10.12) 
(iv) 
(10.13) 



\h£U m {x, 0; y, T) - h~)U m ,(x, 0; y,T)\< ch 2 m 
\h^Vh m U m (x, 0; y, T) - K$V hm , U m > (x, 0; y,T)\< ch 2 m 
\h£b h JJ m (x,tyy,T) - h-)A hm ,U m ,(x,Q;y,T)\< ch 2 m 

KtU m {x, 0; y, T) - h^U^x, 0; y, T)\< ch 2 m 



Proof. It suffices to show this inequality for m' = m + 1. Let us assume for simplicity that 
L = 2 m "h„ l0 and A m — {x — h m i, i — 0, ....2 m } for all m > m . The argument extends to more 
general lattice geometry but the consideration of these more general cases would obscure the 
simplicity of the proof with needless detail and will thus be omitted. 

Let B m be the Brillouin zone as defined in equation (5.161. Let T m : £ 2 (A m ) — > £ 2 (B m ) be 
the Fourier transform operator defined so that: 



(10.14) 



f(p)=T m (f)( P ) = h m E f(x)e 



IpX 



x£A„ 



for all p € B m . The inverse Fourier transform is given by 



(10.15) 



p£B„ 



The Fourier transformed generator is diagonal and is given by the operator of multiplication 



by 

(10.16) 

We have 
(10.17) 



i m {P) = T m L m T m (p,p) = -ifi — - h a 



m 



2L 



E 

p£B m 



hi 



,Tt m (p)Ap{y-x) 



Using this Fourier series representation, we find 
\KnU m (x,0;y,T) - h m 1 +1 U m+ i(x, 0; y, T) 
1 



< 



2L 



E 

p£B m 



j e ip(y-x) 


1 




+ 2L 



E 

p£B m+1 \B„ 



e Te m+1 (p) e ip(y-x) 



(10.18) 
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Let 

(10.19) 



K„ 



log h 



m+l I 



<7 2 T 



If h m is small enough, i.e. if too is sufficiently large, we have that 
1 



2L 



E 

peB m ,\ P \>K„ 
(10.20) 



e Tl m (p) e ip(y-x) 



2L 



pefl m ,|p|>x„ 



where 3?(a) denotes the real part of a G C and c denotes a generic constant. Similarly 

1 



1 

2L 



E 

P eB m+1 ,\ P \>K„ 



e T£ m+1 (p) e ip(y-x) 



< _ V e TK(l m + l {p)) 

2L ^ 



< cexp Tu 



peB m ,\ P \>K m 
2 cos h m+ iK — 1 
'V+i 



< ch 2 m+1 



(10.21) 
Since 

(10.22) 

and 

(10.23) 



, , .i i sin hp sin2hp 1, 9 q 
h 2 p 3 - -h 4 p 5 < < -h 2 p 3 



1 2 4 COS hp -1 COS2/ip-l 1 24 1 4 6 

8* P (2h)^--8 hp + 48 kp - 



we find that if \p\< ^ then 
(10.24) 

Moreover, since 
(10.25) 



\L(p)-L+i(p)\< jh 2 \p\*+^hY- 



1 9 COS hp — 1 In 1,94 

<--p 2 + -fey 



-jr < 

2^ - h 



2" 24 



we conclude that in case \p\< h the following inequality holds: 



(10.26) 

Hence, if toq is large enough, we find 



cos hp — 1 In 

< — p 

h ~ V 



1 

2L 



E 

pes ro ,| P |<if 



,Tt m (p) _ TZ m+1 (p) \ p ip(y-x) 



(10.27) 



^ peB m ,\ P \<K 

- E 



2L 



e 4 



pGB m ,\p\<K 



p2 (^-hl\p\ 3 +^h 2 m p>) -irhi 



< 



for some constant c > independent of to. This concludes the proof of the bound in ( 10.10 1. 
To estimate the sensitivity in (10.11) notice that 

(10.28) 



h" 1 Vl7 m (x,0;tf,T) = 7 E e 



T£ m (p) smph m ip(y- x ) 



L 



peB„ 
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and 



\Kn U m(x,0;y,T) - h rt L 1 U m+ i(x,0;y,T)\ 



< 



2L 
1 

2L 



E >< 

p£B„ 

E 

P eB m+1 \B„ 



TZ m (p) 



sin ph 



i(p) smph m+1 ^ rlp(y _ x) 

h-m+l 



e Tl m+1 (p) e ip(y-x) 



(10.29) 
Let 

(10.30) 



K„ 



| log h 



m+l\ 



a 2 T 



If h m is small enough, we have that 



1 

2L 



E 

peB m ,\ P \>K r 



c Ti m (p) Smpfem c ip(y-x) 



< c 



< 



2L 



E 



^T5R(i m (p)) sm P^m 



P6-B m ,|p|>-KT„ 



sin if /i„ 

fom. 



exp TV 



, cos i<T/i m — 1 



(10.31) 

where c denotes a generic constant. Similarly 



(10.32) 



1 

2L 



E 

p&B m+1 ,\p\>K T , 



e Te m+1 (p) e ip(y-x) 



< ch 2 . 



If m is large enough, we also find 

' sin ph 



1 

2L 



E 

pes m ,| P |<K 
1 



™ e Tt m (p) 



sin ph n 



hm+l 



■+l p Te m+1 (p) \jp(y-x) 



< 



2L 



E 



peB m ,\ P \<K„ 



hm+l 



sinph m 


e 


ip 2 ^hi l \ P \ 3 +4fh 2 m p i 


h m 






sin ph 


m 


< ch 2 m 


h m 







(10.33) 



for some constant c > independent of m. This concludes the proof of the bound in (10.11 1. 
The bound in (10.12) can be derived in a similar way. 

Finally, consider the following Fourier representation for the discretized kernel 

(10.34) h m l U^{x^;y,T) = ± £ (l + 

Consider the formula 



peB n 



(10.35) 



l + <5*4(p) =exp(riog(l+4(») 
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and let's represent the difference between kernels in (10.131 as follows: 

\h- x U m (x, 0; y, T) - h^U%(x, 0; y, T) \ 



< 



2L 



peB„ 



( exp (T£ m (p)) - exp ( § log (l + 6U m (p)) 



St 



,ip(y-x) 



< 



2L 



L E 

eB m , P <K m 

- E 



T 



exp ( - log (1 + 8U m (p)) - T£ m {p) ) - 1 



p£B m ,p>K„ 



exp (T£ m (p)) 



- E 

2L ^ 



p£B m ,p>K„ 



T 



exp ( — log (1 + 8U m (p)) 



(10.36) 



where K m is chosen as in (10.19). The very same bounds above lead to the conclusion that this 
difference is < ch^. 

□ 



11. Kernel Convergence Estimates for Diffusion Processes 

Diffusions are a particularly important class of Markov processes which generalize Brown- 
ian motions to allow for space dependent drifts and volatility. In this section we find kernel 
convergence estimates for the one-dimensional case following (Albanese 20076). 

Let A = (A m ), m — mo, mo+1, ... be a simplicial sequence converging to an interval linim^oo A„ 
[-L, L] C R, where < L < oo. Let /i(x) and tr(x) be smooth functions defined in a neighbor- 
hood of [-L, L]. The generator of a diffusion on A m has the form 

(11-1) An = A»(a;)V hm + ^a(x) 2 A hm . 

We assume that mo is large enough so that 

a 2 (x) ^ | M (x)| 



(11.2) 



2hl 



> 



2h. n 



for all m > mo and all x £ A m . The definition of the generators at boundary points can be 
extended by imposing one of the boundary conditions in the previous section, i.e. reflecting, 
absorbing, periodic or mixed. 

Theorem 17. (Convergence Estimates for Diffusions.) Consider a diffusion process de- 
fined as in {ll.ty where fi(y) and o~(y) are smooth functions satisfying the bounds in (11.2) for 
all m > mo- Assume that boundary conditions are either periodic or absorbing. Then there is a 
constant c > such that 

(11-3) \hnU m {x, 0; y, T) - h m )lJ m ,{x, 0; y,T)\< ch 2 m 

for all m! > m and all y G A rn . 

It suffices to establish the above inequality in the case ml = m + 1. In fact, given this 
particular case, the general statement can be derived with an iterative argument. Let h — h m+ i 
so that h m = 2h. 

Recall that the path integral representation for a diffusion process has the form in equation 



pl.4p. In the special case of a time-homogeneous process, this expansion reads as follows: 



(11.4) 



E 



U m (x,0;y,T)=J2^~ q 

9=1 7 € r m : 7o = x,7 9 = y 
l7j — 7i-il= 1 Vj > 1 
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Re(o) 



Figure 1. Contour of integration for the integral in (11.91. C+ is the countour 



joining the point D to the points E, A, B. C_ is the countour joining the point 
B to C to D. 



where 
(11.5) 



W m ( 1 ,q,T)= [ ds x [ ds 2 ... [ ds q e^- s ^ c ^y^l[ fe^+ 1 - s ^ £ -^^)£ m ( 7i) 7,-+i) 

JO J 81 J Sg-l j_Q \ 

where sq = 0. 

Let us introduce the following two constants characterizing the volatility function: 



(11.6) 

and let 
(11.7) 



E = inf <r(x), Si = sup \J <j(x) 2 + h m \fi(x)\. 

xeA m xeA m 



M = sup \n(x)\. 

x£A m 



Since our interval is bounded, we have that Eo > and Si,M < oo 
Let us introduce the following Green's function: 

(11.-' ' ' 



G m (x,y;uj)= / U m {x, 0; y, t)e luJt dt= . (x,y). 



The propagator can be expressed as the following contour integral 
(11.9) 



U m (x,0;V,T) = [ ^G m (x,y;co)e^ T + f ^G m (x,y;Lo)e^ T . 
JC- 2lT Jc+ 27r 

Here, C+ is the contour joining the point D to the points E, A, B in Fig. [T| while C_ is the 
contour joining the point B to C to D. By design, each point to on the upper path C+ is 
separated from the spectrum of C. 
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Lemma 1. For m sufficiently large, there is a constant c > such that 

duo 



(11.10) 



2tt 



G m (x,y;uj)e 



< ch 6 



Proof. The proof is based on the geometric series expansion 

1 ^ 1 



(11.11) 



G m (u>) = 



3=0 2 



he 2 A" 



whose convergence for ui £ C + can be established by means of a Kato-Rellich relative bound, see 
(Kato 1976). More precisely, for any a > 0, one can find a (3 > such that the operators V m 
and A m satisfy the following relative bound estimate: 



(11.12) 



| 2 <a||A"7|| 2 +/3||/||, 



for all periodic functions / and all m > rriQ. This bound can be derived by observing that V m 
and A m can be diagonalized simultaneously by a Fourier transform, as done in the previous 
section, and by observing that for any a > 0, one can find a [3 > such that 



(11.13) 



sin h m p 



< a 



cos h m p — 1 



hi 



for all m > mo and all p € B m . 

Under the same conditions, we also have that 



(11.14) 

Hence 

(11.15) 



\^ m f\ 



2Ma 



ct 2 A"7 



— / 



< 



2Ma 



V 2 



— <r 2 A m -j - 

2 i(7 2 A m 



PWSh 



-f 



■p 



ia 2 A r ' 



— f 



< 1 



where the last inequality holds if lu G C+, if a is chosen sufficiently small and if m is large enough. 
In this case, the geometric series expansion converges in (11.11 1 converges in I? operator norm. 
The uniform norm of the kernel \G m (x,y;w)\ is pointwise bounded from above by ft." 1 . 

Since the points B and D have imaginary part equal at height 4 l lo s^ lm l ; the integral over the 
contour C + converges also and is bounded from above by ch^ in uniform norm. 



□ 



Lemma 2. If q > e 2 ^ 2 lT we have that 



(11.16) W m (%q;0,T) < 

Proof. Let us define the function 



exp 



£ 2 T 



(11.17) 



T 2 



2h 2 



l(t > 0) 



where \(t > 0) is the characteristic function of R+. We have that 
(11.18) W m ( % q;0,T) <4>* q {T) 

where (j)* q is the q— th convolution power, i.e. the q— fold convolution product of the function 
by itself. The Fourier transform of <fi(t) is given by 



(11.19) 



2/ii 



2iwh 2 



The convolution power is given by the following inverse Fourier transform: 



(11.20) 



(t>* q (T)= / (j)(w) q e iulT = 



Si 
So 



2</ 



2iw/j 2 , 

V2 
^0 



2?r' 
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Introducing the new variable z = 1 



2iuihi 



the integral can be recast as follows 



(11.21) 



(j)* q (T) = 



^2-2(j v 2<? 



lim 



■ exp 



2/a 



(z - 1) d 



where Cr is the contour in Fig. [2] Using the residue theorem and noticing that the only pole of 
the integrand is at z = 0, we find 



(11.22) 



S 2 T 
(9-1)! 



exp 



Making use of Stirling's formula ql 

q 



2irq q+ ze q , we find 



(11.23) 



If log q > log 



(f>* q (T) 



exp 



EgT 
2/ii 



Q log 



-S§T 
2hl 



+ q(l-logq) 



2, then we arrive at the bound in (11.161 



□ 



Definition 41. (Decorating Paths.) Let m > mo and let 7 = {yo, y\, yi, •■•■} be a symbolic 
sequence in T m . A decorating path around 7 is defined as a symbolic sequence 7' = {yo, y'1,1/21 ■••■} 
with y[ G h m+ iL containing the sequence 7 as a subset and such that if y'j — yi and y' k = yi+i, 
then all elements y' n with j < n < k are such that \y' n — y'A< Kn+i- Let V m+ i(^) be the set of 
all decorating sequences around 7. The decorated weights are defined as follows: 



(11.24) 



W m (7,3;0,T) 



OO 

E £ 



W m+l {i,q';Q,T). 
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Finally, let us introduce also the following Fourier transform: 



(11.25) W m ^,q;cu) = / VK m ( 7 ,g;0,i)e 1 ^, W m ( 1 ,q;uj)= / W m ( 7 , q; 0, t)e Mt dt. 

Jo Jo 

Definition 42. (Notations.) In the following, we set h = h m+ i so that h m — 2h. We also 
use the Landau notation 0(h n ) to indicate a function f(h) such that h~ n f(h) is bounded in a 
neighborhood o/(0). 

Lemma 3. Let x,y € A m and let C_ be an integration contour as in Fig. [7] Then 

n dw 



(11.26) 



2G m+ i(x, y; lu) - G m (x, y;u>) e 



2tt 



= 0{h A ). 



Proof. We have that 
(11.27) 

OO 

2G m+ i(x,y;cj) - G m (x,y;uj) = ^2~ 9 ^ (2W m (j , q; lu) ~ W m (-f , q; lu) 

9=1 7 G r m : 7o = x, j q = y 
hj -7j-i|= 1V J > 1 

The number of paths over which the summation is extended is 



(11.28) N(~/,q;x,y) = ${l G r m : 7o - x, lq = y, | 7j - 1] _ 1 \= lVj > 1} 

where fc = ^ . Applying Stirling's formula we find 

(11.29) 
Hence 

2G m+1 (x, y;uj) - G m (x,y;oj) )e lujT 



q 

l + k 



N-y < 2 q 



nq 
2^ 



< c ti - max 

q=1 V q 7 g r m : 7o = x, 7g = y 



2W m ( 7 ,g;w) - # m ( 7 ,q;w) e 



2^ 



It? -7i-il= 1v j > 1 



(11.30) 

for some constant c « ^/|~ > 0. It suffices to extend the summation over q only up to 

_ e 2 E?T 



(11.31) 



2h 2 



To resum beyond this threshold, one can use the previous lemma. More precisely, we have that 

i dbJ 



2G m+ i(x, y; lu) - G m (x, y;u>) e 



2tt 



2W m (7,g; w) - W m ( 7 ,«?; e 



i da; 



2?r 



< <V<7max max 

9J6 T m : 7o = a;, 7? = y 

iTj -7?-i|= lVj > 1 

(11.32) 

Let v(x) = a{x) 2 . To evaluate the resummed weight function, let us form the matrix 



/ 



v(x-\-h) 



v(x-\-h) fi(x-\-h) 



n 
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and decompose it as follows: 

(11.34) £(x; ft) = pA)(x) + \^ x ) + £ 2 (x) + h£ 3 (x) + 0(h 2 ). 
where 

(-v(x) \v{x) \ 

(11.35) C (x) = \v(x) -v\x) \v{x) , 

V \v{x) -v(x)J 

-v'(x) \v'(x) - \»{x) 

(11.36) £i(x) = \ \n{x) -\»{x) 

■\v'{x) + \ 



-\v'{x) + \n{x) v'(x) 



-±v»(x) iv"(x)-fr'(x) 
(11.37) C 2 (x) = ( 

\v"(x)-y(x) -\v"{x); 

and 



6 



v"'{x) ±v"'(x) - \n"{x) 



(11.38) C 3 (x) = I 

-^ {x ) + y'{x) \v'"(x) / 

Let us introduce the sign variable r = ±1, the functions 

(11.39) <fo{t,x,T) = 2£ m (x,x + 2Th)e tc ^ x < x h(t > 0) 

(11.40) <l>i{t,x,T) = 2C m+1 (x + Th 1 x + 2Th)e t£(x ' h) (x 1 x + Th)l(t > 0) 
and their Fourier transforms 
1 , x ( v ( x ) K x )\ ( v ( x ) 

K», *. r) = + + V " iX)+ / iX) + {^ + ^)rk + 0(tf ) 

(11.41) 

< x| (— £(x; ft) + iw) |x + rft > . 

where 

(11.42) |x>= 1 , and |x + rft>= 

W \<y T ,-i y 

We also require the functions 

(11.43) V>o(M) = e t£ "*( x ' x h(i>0), ViM ^e^^^l^O) 
and the corresponding Fourier transforms 

(11.44) V>o(w,x) = ( +iu)j , Vu(w,x) =< x\(-C(x; ft) + iw) |x > 
If 7 is a symbolic sequence, then 

9-1 ^ 

(11-45) Wm(7,9;w) = Mw,7 9 ) JJ ^ (w;7j,sgn(7 j+ i - 7,)) 

9-1 ^ 

(H-46) T^ m (7,g;w) = $i(w,7,) JJ ^i(w; 7^, sgn(7 j+ i - 7,-)). 
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Let us estimate the difference between the functions <pi(u>, x, r) and 4>2(^, x, r) assuming that 
u> is in the contour C_ in Fig. [2J Retaining only terms up to order up to 0(h 3 ), we find 

2 , n 1 2fi(x)rh Aiujh 2 o . .iurh 3 16w 2 /i 4 ^,, 5 , 
0o(w, s,r) = 1 + ry - y^ - —r^- + 0(h 5 ). 

(11.47) 

A lengthy but straightforward calculation which is best carried out using a symbolic manipulation 
program, gives 



i(uj,x,t) = 1 + 



2fj,(x)rh 4iajh 2 



v(x) v(x) 



[8fM{x) - v'(x)] 



rh 3 



v(x) 2 

1 A, .2 l4 

+ r{x) ■ h 3 r + iujh 4 p(x) %r + °( h5 ) 

v{x) z 



(11.48) 
where 

r(x) 

p(x) 
(11.49) 



1 



2v{x) 3 
1 



[n"(x)v(x) - 4[i(x) 3 + 2v'(x)fi{x) 2 - 2v'(x)v{x)fi'(x) 

~ {2n{x)ix'{x) + v"{x)v{x) - 2v'(x) 2 ) f i(x)] . 



[4/z(x) 2 - 2v'(x)n(x) + Av(x)n'{x) + v"(x)v{x) - 2v'(x) 21 



v(x) s 



We have that 



( log^o(w;7i)Sgn(7i+i - 7i)) ~ l°g<£i( w ;7i> sgnf^+i - 7,-)) j 



q-l 



E f^|^ + r ^)) ^frj+i ~ 7i) + (|o;||b|| 00 +2|a;| 2 ||«- 2 || 00 )0(/ l 4 g) 



«(7g) 



^(7o) 



= iw/1 2 log 

(11.50) 

where i?(x) is a primitive of r(x), i.e. 
(11.51) R(x) = I r(z)dz 



+ h 2 (R( lq ) - i?( 7o )) + (|a;|||p|| 0O +2| w | 2 || W - 2 || oo )0(/i 4 g) 



We conclude that there is a constant c > such that 

.q-\ q-l 

(11.52) 



/ ( II ^o(^;7i,sgn( 7j - +1 - t,-)) - Y[ 0i(w,7j,sgu(7 i+ i - 7j )) )e l " T 



dw 
2^ 



< c/1 2 



for all q < q max . Here we use the decay of e lulT in the upper half of the complex oj plane to offset 
the u> dependencies in the integrand. Similar calculations lead to the following expansions: 

Ah 2 2h 2 1 

(11.53) M",y) = tt + oK), M",v) = tt + oM 4 ) = -M^,y) + o(uh 4 ). 

v{y) v(y) 2 



Since q < ch 2 and lo < |log/i|, we find 



(11.54) 



2G m+1 (x,y;uj) - G. m (x,y;uj) e 



AuiT 



2^ 



< cq m ^h 4 < ch 3 . 



This completes the proof of the Lemma and of the Theorem. 



□ 
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12. Convergence of Time Discretization Schemes 

In this section we analyze the convergence of the time discretized kernel that is obtained by 
means of fast exponentiation. 

Let A — (A m ), to = too, mo+1, ... be a simplicial sequence converging to an interval lirn m _,oo A m — 
[0, L] C R, where < L < oo. Let fi(x) and cr(x) be functions on [0, L] satisfying all the condi- 
tions in the previous section and consider the generator of a diffusion on A m of the form 

(12.1) C m = Li{x)V hm + -a(x) 2 A hm . 

Theorem 18. (Convergence Estimates for Fast Exponentiation.) Let St > and con- 
sider the discretized kernel 

(12.2) U%(x, 0; y, T) = (1 + 6t£ m ) % (x, 0; y, T). 



where C m is the operator in (12.1) and St m is so small that 

(12.3) min 1 + 5t m £ m (x, x) >0 

xeA m 

Assume that boundary conditions are either periodic or absorbing and that the ratio £ = N is 
an integer. Then there is a constant c > such that 

(12.4) \K^U m {x,Q;y,T) - h^U%(x,0;y,T)\< ch 2 m 
for all to > toq and all y € A m . 



Proof. A Dyson expansion can also be obtained for the time-discretized kernel and has the form 

oo N N N 

i&w^rHE E EE- E 

q=l 'Y£T m :~fo=x,~f q =y fci=l fe 2 =fci+l fc g =fc 9 -i+l 

(12.5) (l + aC™(7o,7o)) 1 (<5t)«[] /: ™(7 J -i>7,)(l + ^ m (7 J ,7,) X ^ 

where t q+ i = T and k q +\ — N . In this case, the propagator can be expressed through a Fourier 
integral as follows: 

(12-6) 0; 2/ 2 , T) = T G«(y 1)2/2 ;a,)e^ 

J- ft 7r 

where 

T_ 
St 

(12-7) G 5 „ t l (y 1 ,y 2 ;cj) = St^U^A V2, jSt) e -^ 5f . 

3=0 

The propagator can also be represented as the limit 



(12-8) U£( yi ,0;y 2 ,T)= lim / G%{y u y 2 ;u)e 

H ^°°JCh 27T 

where C// is the contour in Fig. |3] This is due to the fact that the integral along the segments BC 
and DA are the negative of each other, while the integral over CD tends to zero exponentially 
fast as 9(w) — > oo, where 3(w) is the imaginary part of u. Using Cauchy's theorem, the contour 
in Fig. [3] can be deformed into the contour in Fig. [T] To estimate the discrepancy between the 
time-discretized kernel and the continuous time one, one can thus compare the Green's function 
along such contour. Again, the only arc that requires detailed attention is the arc BCD, as the 
integral over rest of the contour of integration can be bounded from above as in the previous 
section. 

Let h = h m and let us introduce the two functions 

(12.9) <M*,a;,T) = 2C rn (x,x + Th)e tCm{x - x) l(t > 0), 

(12.10) cj> st (j, x, t) = 2C m (x, x + rh) (1 + 5tC m {x, x))^ 1 . 
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\E-%/St 
D 



lm(co) 



-7T/6t 



c/h 



iH+w/6t 
C 



Tf/St Re ( f0 ' 



Figure 3. Contour of integration for the integral in (12.81 



and the corresponding Fourier transforms 

--s i \ f°° , , \ 4,,fduj ( v(x) fi<(x)\ ( v(x) 

(12.11) <fa(w,x,T)= Mt,x,T)e~ 1 ^ -' 



2tt 



h 2 



h 2 



(12.12) fe(u,i,T) 



We have that 



5t(w,x,r) 



i/(x) ( _fJ,(x) 



u(x) oj 2 



<~ — ; — I I iw + — tV" — -r-ft + <3(<5£ 2 ) 



(12.13) 



0o(w, x, r) + -^h 2 St + 0{h 2 8t 2 ). = o (w, x, r) + (9(/i 4 ), 
2v(x) 



where the last step uses the fact that St — 0(h 2 ). 
Let us also introduce the functions 



(12.14) <Mt,x,r) = e t£ "^l(t>0), ^t(j, x, r) = £ (l + StC m (x, x)) 



3-1 



fc=l 



and the corresponding Fourier transforms 

-l 



(12.15) ifotw.z.T) 

Again we find that 
(12.16) 



v(x) 



10J 



ii) St {w,x,T) = ( e iuSt -l + 5t 



h 2 



4>q{u,x,t) = 4> St (LJ,x,r) + 0(h 4 ). 
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If 7 is a symbolic sequence, then let us set 



(12.17) 
(12.18) 

We have that 
(12.19) 



9-1 ^ 

W m {l, q; = -0o(w, 7,) J| o (w; 7j,sgn(7 i+ i - 7,-)) 
j=o 
9-1 ^ 

W^*(7,g;w) = ^5t(w,7g) II 05t(w;7j,sgn(7 i+ i - 7j-)). 



G%(x,y\w) - G m (x,y;u) = £V 



E 



9=1 



W*(7,«;w)-^ m (7,g;w)) 



7 G r m : 7o = 2;, 79 = V 
|Tj -7j-i|= IVj > 1 

The integration over the contour in Fig. [T] can again be split into an integration over the 
countour C_ and an integration over C+. The integral over C+ can be bounded from above 
thanks to Lemma [3] Furthermore, we have that 

1 duj 



G S ^{x,y;uj) - G m (x,y;uj) ]e 



2tt 



(12.20) 



< <V<j max max 

9, 7 G T m : 70 = x, 7, = y 

\lj -lj-i\= IVj > 1 

< c/i 3 



W^( 7 , ? ;a;)-W m (7,g;a;) )e 



2n 



□ 



13. Hypergeometric Brownian Motion 

This section is based on work in collaboration with Joe Campolieti, Peter Carr and Alex 
Lipton, see (Albanese et al. 2001). 

In this section we pose the problem of constructing driftlcss diffusion models which reduce to a 
given diffusion by means of a combination of a measure change and a coordinate transformation. 

Consider a Markov process on the simplicial sequence A m C R d with generator C m (x,x'\t). 
Let p be a real valued parameter and suppose fm( x ) and fm( x ) are two linearly independent 
solutions of the equation 



(13.1) 



^2 £rn(x, x'; t)f m )(x') = pf m (x) 



x'£A,- 



for all x 6 Int(A m ). Notice that on the boundary d(A m ) these equations may fail and, in actual 
applications, they will indeed as a rule fail. Consider the function 

(13.2) g m (x; t) = e^d/^) + c 2 f m {x)) 

for some choice of constants c\,C2 such that this function is strictly positive. This function 
satisfies the equation 

dg m {x;t) 



(13.3) 



dt 



+ (C m g m )(x;t) = 0. 



Hence g m (x;t) defines a measure change and one can construct a new Markov generator by 
setting 

1 dg m (x;t) 



(13.4) 
Consic 1 
(13.5) 



rg ( t. j.\ 9m{% j t) p , / ,\ 

Z-, m yXj X , L) — - — L, ra \X^X , Z) ~r 



9rn(x;t) 

Consider the linear fractional transformation 



g m (x;t) dt 



(^) 



Clfm( x ) + °2fm( x ) 
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for some choice of constants 03,04. 



Theorem 19. (Linear Fractional Transformations.) The process Y(xt) satisfies the mar- 
tingale condition 



(13.6) 

for all Xt G Int(A m ). 
Proof. We have that 



lim -E t 

sj.0 S 



Y m (xt+ S ) — Y m (xt) 



= 0. 



lim s i \E t 



Y m {x t + B )-Y m (x t ) 



^x'£A m £m( X ' X '){Ym{x') Y m (x)) 



^T)(j:^A m C m (x,x')(c z fl(x') + c 4 /^'))) + g J x t)2 d ^ (c,fl(x) + cj« (a)) 



- at 



e- pt (c 3 /, 1 „(^)+C4/^( a; )) 



= 



(13.7) 



□ 



This Theorem provides a general methodology for constructing a process which is nearly a 
martingale out of a Markov process. In the particular case of one-dimensional diffusion processes, 
the construction gives rise to families with up to 7 adjustable parameters of analytically solvable 
diffusions with drift equal to zero within the interior of the domain of definition. What makes 
the case of one-dimensional diffusions special is the fact that the function Y m (x) is invertiblc in 
the limit as m — > 00, i.e. either monotone increasing or monotone decreasing in this limit. 

More specifically, consider the diffusion with Markov generator: 

(13-8) C m = 8{x) V hm + ^(x) 2 A hm . 



where i£E. Two important cases are the CIR process (reducing to the Bessel equation in the 
continuous limit) for which 

(13.9) 5(x) = (Ao + Xix), v(x) — v^\fx 

with x € K+ and the Jacobi process (reducing to a gaussian hypergeometric polynamials of type 
2-Fi) for which 

(13.10) S(x) = (Ao + Aix), v(x) = VQ^/xJl - x) 
and x G [0, 1]. 

The construction above admits a continuous limit if the simplicial functions /^(x), j — 1, 2, 3, 4 
converge to twice diffcrcntiable functions P(x),j = 1,2,3,4 and satisfy the equation 

(13.11) Cf = pf 
within the interior of the domain A = linim—.oo A m . Here 

d 



(13.12) 



C = 5(x) 



dx 



r {x) a*- 



Theorem 20. (Invertibility.) Let f J (x),j = 1,2,3,4 be functions satisfying equation 13.12 
the interior of the corresponding domain of definition and let 

c^f 1 {x) + c i f 2 {x) 



(13.13) 



Y{x) 



cif 1 (x) + c 2 p(x)- 

for some choice of constants Ci, C2, C3, C4 such that the denominator in this equation has no zeros. 
Then we have that 

W{y) 



(13.14) 



Y(x) = 



g(y)' 



7 dy + const 
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where 
(13.15) 



W(x) = ±cr exp ( - J 



x 2S(y) 
v{y) 2 



dy 



where either sign would be allowed and o~o is a positive constant. In particular, the function Y(x) 
is invertible. 



Proof. Let us introduce the Wronskian 

(13.16) 

such that 
(13.17) 



„„ , dh(x) , . dg(x) , . . 



dx 



dx 



g{x) = cif{x) + c 2 f 2 (x), h(x) = c 3 f{x) + Ci f 2 (x). 
A direct calculation shows that 

d 



(13.18) 
and 
(13.19) 



dx 



Y(x) = W(x)Y(x) 



dx 



W(x) 



v\xy 



Hence, ( 13.15| ) gives the solution to the equation in (13.191. 
Let X(y) be the inverse of the function Y(x), we have that 

(13.20) 

where 



□ 



dy =± _dX(y) 



(13.21) 



a(y) = Y'(X(y)HX(y)) 



°{y) HX(y)) 



g{x{y),p) 2 

Theorem 21. (Kernel mapping.) In the limit m — > oo, the propagator density for the process 
Y t with zero drift and volatility as in (13.21) is given by 

(13.22) 



U(y , 0; y, t) = 4^A^(X(y )A X(y),t) 



<r(y) g( x (yo),p) 

where u(xq,0; x,t) is the propagator density for the process of generator fcl3.12\ i. 



Let us consider in detail the case of the CIR model in (|13.9|). If Ai = 0, then the pricing 

functioi 

Ay/XXQ " 



kernel for the state variable is expressed in terms of modified Besscl functions as follows 

hi^r- 1 ) p -2(x+x )/v$t 





(13.23) 



u(x, i; xq, 0) 



x 



Xq 



i#/2 



I2 



vlt 



The generating function is 



(13.24) 



v(x, p) = x 



' 8px 



/ 8px 



with arbitrary constants gi,<?2- Here I v (z) is the modified Bessel function of order v and K v {z) 
is the associated McDonalds function. In this case we obtain a dual family with 6 adjustable 
parameters. 

In case Ai < 0. the propagator density for the state variable x can still be expressed in terms 
of modified Bessel functions as follows: 



(13.25) u{x, t; xq, 0) = Ct 



>-Ait\ K^-l) 



:r 



exp 



[-c t (x e Xlt + x)] l2x SL _ 1 ( 2ctVxx e x ^ 
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where ct = 2Ai/(z^o(e Alt — 1)). For a derivation see (Giorno et al. 1988). The general solution of 
equation ( 13.1 ) reduces to Whittaker's equation and generating functions have the general form 

2Ai 



(13.26) 



v(x,p) 



q{W k , 



2Ai 



o / v "0 

for arbitrary constants 91,92- Here Wk, m (') and Mfc im (-) are Whittaker functions which can also 

be expressed in terms of confluent hypergeometric functions or in terms of Kummer functions. (Abramowitz 

and Stegun 1972) This construction gives rise to a dual family with 7 free parameters where 

Ao , p A n 1 



(13.27) 



k = 



Ai 



Ao 
, ,2 



The 7 parameter family which reduces to the CIR model has a local volatility function defined 
on either an interval or on a half line and behaves asymptotically as the CEV volatility on one 
hand and as a quadratic model on the other. This hybrid shape allows for greater flexibility. 

Next, we show that classic exact solutions in the literature, namely quadratic and CEV 
models, can all be rediscovered as particular cases of our general formula for the Bessel family 
where we make use of the above solutions to the underlying x space process with (3 = §, Ai = Q 



and A = Ao. 
where 

(13.28) 



Without loss of generality, we can fix v$ = 2. Let's specialize further to the case 



y 4-i( 




which leads to a transformed process yt — Y(xt) with volatility 
(13.29) a(y) = 9 



where x — X(y) is the inverse of the function in equation 
are positive, y is arbitrary and A > 2 



113.28 1 



In this family, a and p 
The function Y(x) maps the half line x € [0,oo) 
into y G (— 00, y], where Y{x) is a strictly monotonically increasing function with dY(x)/dx = 
<j{Y(x))/v(x). This solution region can be inverted so that y € [j/,00). This is accomplished 
by either replacing a by —a in equation (13.281 or by applying a linear change of variables that 
maps y into 2y — y. In this special case, we make use of the generating function in equation 
(13.241, with the choice 92 = 0, and formula ( 13.22| reduces to 

<i 3 . 3 o) [ , Mfcl) .L^^««: i , (vmms), 

at J | _ 1 (V2pX(y )) 5 \ t J 

We note that this density integrates exactly to unity in y space (i.e. no absorption). 

Example 4. (Quadratic volatility models) Pricing kernels for quadratic volatility models 
are readily obtained as a subset of the above general family with the special choice of parameter 
A = 3. After making the substitution y — ► 2y — y and setting a — ( y — y)/n the transformation 
function Y{x) becomes 



(13.31) 



rW-» + »- B ** ( "rff -' + -r5& 

7T ll(<Toy/x/2) exp((70V a; ) 



where do > 0- Here, we assume that y > y. The inverse transformation X(y) is given by 
(13-32) X(y) = (l/a 2 ) log 2 [l + (y - fj)/(y ~ y)}, 



Bessel function of order | 



and the volatility function a(y) is obtained by insertion into equation (13.291 while using the 
Bessel f 

(13.33) 



00 



= ""=, (y -y){y- y)- 
{y-y) 
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Inserting the expression (13.321 into equation (13.301, one obtains the pricing kernel 



(13.34) U( y> t;y o ,0) = J {ya ^ ^^^^0?)/^^ ( 



a(y)V2nt)l (y - y)(y - y) "V °o* / 

where <j>{y) = log((y — y)/(y — y)). In the special case of a volatility function with a double root, 
i.e. 

(13.35) a(y) = a a {y - y) 2 

the pricing kernel is computed by taking the limit as y — > y, and one finds 



U(y,t;y o ,0) 



CT ' 



J (yo - y) 

faftiy-y) 3 



(13.36) 



e -((tf-S) _1 -C»o-S) _1 ) /2<rgt 
e -((l/-fi)~ 1 +(vo-i/)- 1 ) 3 /2(rgi 



Example 5. (Lognormal models) The pricing kernel for the log-normal Black-Scholes model 
with a{y) = a$y is a particular case of the above formula for the quadratic model. The derivative 
with respect to y of the quadratic volatility function in (13.331, evaluated at y — y, is <tq. 
Taking the limit y — > —00 (or y << y), while holding the other parameters fixed, one obtains 
a(y) = (T{){y — y). The pricing kernel in ( 13.341 gives the kernel for the log-normal model in the 
limit y — > —00, i.e. 

2 / 



(13.37) U(y,t;y ,0) = 



1 



(y - y)oo\/27rt 



exp 



-tlog((y -y)/(y-y))-^-t 



2ah 



Example 6. (CEV model) The constant-elasticity-of-variance (CEV) model is recovered in 
the limiting case as p — > 0. Assume A > 2 and let 9 > be defined so that A = 0~ l + 2. The 
transformation y — Y(x) 

(13.38) Y(x)=y + (a 2 x)-W- 1 
has inverse x = X(y) given by 

(13.39) X(y)=a^(y-y)- 2e , 
for any constant y. The volatility function for this model is 



(13.40) 



°{y) 



(y-y) 



1+0 



In the limit p — > 0, the Laplace transform v(X(y),0) — 1, which implies that the numeraire 
change is trivial in this case. The pricing kernel can be evaluated by substitution into the 



general formula (13.221, and after collecting terms, it turns out to be 

\0\ (yo-V)5 r -((y-y)- 2e +(y -y)- 2e )/2a 2 n t 



U(y,t;yo,0) 



alt (y 



y) 



\+28 



(13.41) 



I 



((y - y)(yo - y)) 



a 2 t 



This formula was derived in the case 9 > 0, for which the limiting value y — y is not attained 
and the density is easily shown to integrate to unity (i.e. no absorption occurs and the density 
also vanishes at the endpoint y — y). We note that the same formula solves the propagator 
equation for 9 < 0, leading to the same Bessel equation of order ±(26>) _1 . In the range 9 < 0, 
however, the properties of the above pricing kernel arc generally more subtle. In particular, one 
can show that the density integrates to unity for all values 9 < —1/2, hence no absorption occurs 
for 9 G (—00, —1/2). The boundary conditions for the density can be shown to be vanishing at 
y — y (i.e. paths do not attain the lower endpoint) for all 9 < —1. In contrast, for 9 S (—1, —1/2) 
the density becomes singular at the lower endpoint y = y (hence this corresponds to the case 
that the density has an integrable singularity for which paths can also attain the lower endpoint, 
but are not absorbed). For the special case of 9 = —1/2 the formula gives rise to absorption. 
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[Note that only for the range 9 e (—1/2,0) the above pricing kernel is not useful since it gives 
rise to a density that has a non-integrablc singularity at y — y. In this case, however, another 
solution that is integrable is obtained by only replacing the order (2#) _1 by — (20) _1 in the 
Bessel function. The latter solution for the density does not integrate to unity and hence gives 
rise to absorption which can be of use to price options in a credit setting.] The special case of 
9 = — 1 gives a nonzero constant value at the lower endpoint, and recovers the Wiener process 
with reflection and no absorption on the interval [y, oo) with 

(13.42) U(y, t; y , 0) = X -= ( e -(»-»>) a /2<xS* + e -(v+yo-2vf/2^t 

14. Stochastic Integrals for Diffusion Processes 

This section provides a new derivation of a theorem by Cameron, Feynman, Girsanov, Ito, 
Kac and Martin, see (Cameron and Martin 1949), (Ito 1949), (Feynman 1948) and (Kac 1948). 
This result is at the center of stochastic calculus and in the following sections, we derive far 
reaching extensions and applications. 

Theorem 22. (Cameron-Feynman-Girsanov-Ito-Kac-Martin.) Consider a diffusion pro- 
cess on the simplicial sequence A m C K and of generator 

(14-1) C m = n(x,t)V hm + ?^-A hm . 

Consider also the process given by the integral 

(14.2) I t = / a(x s , s)dx s + b(x s , s)ds 

Jo 

where a(x,t) and b(x,t) are smooth functions in both arguments. Let us introduce also the 
function 6(x,t) — J a(y,t)dy. We have that 

(i) Ito's Lemma in integral form. 

(14.3) I t = 6{x t ,t) - 6{x ,Q) + Jt + 0(h). 
where 

(14.4) J t = J (b(x s , s) - ^<t(x s , s) 2 a'(x s , s) - 4>(x s , s)^jds. 

where (j)(x,t) = §i<j>(x,i) and a! (x,i) — J^a(x,t). 

(ii) Feynman-Kac formula. The characteristic function of J t on the bridge leading from 
x to y is given by 

(14.5) 

Eo[j pJ *5{xt - y)\x = x]=Pexp(J (c(s) + ipb(s) - |<T 2 (s)a' \s) - ipj>(s) + 0(/i)))(a:,j/). 

(iii) Ito's Lemma in differential form. Let <p(x, t) be a smooth function in both arguments. 
Then we have that 

lim£; t [s~ 1 (<f>(x t+s ,t + s) - (j)(x t ,t))\x t = x] 
t~ . „s 96 , , , . 86 , a(x,t) 2 d 2 6 , , , 

(14.6) = -£{x,t) + fJL(x,t)-£(x,t) + -^y^-gjffo*) + °(M 

and 

(14.7) limE t [s- 1 {6{x t+s ,t + s)-6(xt,t)) 2 \x t = x] = a{x, t) 2 (j^{x, t)j +0{h m ). 

For all n > 3 instead we have that 

(14.8) ]imlimE t [ S - 1 {6{x t+s ,t+s)-6{x t ,t)) n \x t = x] =0. 

hlO sj.0 
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Proof. Consider a discretization of the process It — (AI)n t . This can be accomplished by 
introducing a simplicial complex of one more dimension. Eventually in the proof, we shall take 
the limit as AI — > while leaving h > constant. 
Consider the lifted generator 

' a(x) ii(x,t)' 



C(x, n; x , n ; t) = 



(14.9) 



2h 2 + 2h 
a(x,t) 2 fi(x,t) 
2h 2 2h 
cr(x,t) 2 



S x ',x+hS \ ri -n- 



S x 'x-hS\ n' — n + 



a(x, t)h 

AI 
a(x, t)h 
AI 



S x ' x S(n - n') + b (Sn',n+i - S n >„) 



h 2 -"•■■"•-v- ■-/ ■ A/ 

where [a] stands for the nearest integer to a. The partial Fourier transform in the n variable of 
this kernel is 



C p (x,x';t) = ^£(a;,0;a;^r^;^)e- ^A/^ " l 

n 

' cr(x,t) 2 (J>{x,t) 



+ 



2h 2 + 2h 
a(x,t) 2 fi{x, t) 



(14.10) 



2h 2 
a(x,t) 2 



2h 



exp 



cxp i 



ha(x, t) 
AI 
ha(x, t) 
AI 



Alpj 6 x > tX+h 
Alp) 5 x > x -h 



h 2 



S x , x S(n-n')+ b{ ^-(e^ AI -l) 



In the limit as AI — > we find 



C p (x, x'\ t) = — ^ — A h + (fi(x, t) - ipa(x, t) 2 a{x, t)) V h 



ipa(x, t)/j,(x, t) — -p 2 a(x, t) 2 a(x, t) 2 — ipb(x, t) + 0(h). 



(14.11) 

Introducing the function (f>(x,t) defined above, we find that 



(14.12) j*"K x rt£(x,x';t)e- il "K x 'rt = C p ( Xl x' ;t) + ipb{x,t)5 xx , - ^a{x, t) 2 a'(x, t)5 xx , + 0{h). 
This equality can be rearranged as follows: 

e- ip ^£(x,x'-,t)e ip ^ x 'V + ip<j>(x,t)S xx , = 

(14.13) C p (x, x'; t) - ipb(x, t)8 xx , + fa(x, t) 2 a'(x, t)S xx , + ip<j>(x, t)5 xx , + 0(h). 
Hence 

Pexp ^ J* £ p (s)ds^J(x,x') = 

e i P (*(s,o)-*(x',t)) Pexp ( f*ds(C(s) -ip(b(s) - \cj( S ) 2 a'{ S ) - <j>(s) + 0(h)))(x,x'). 



(14.14) 

The joint kernel is thus given by 



Pexp^ £dsL(s)J(x,I;x',r) = 
I ^eW-^^'V+^WPexp ( /„' ds(C(s) - ip(b(s) - \a(s) 2 a'(s) - 0(a) + 0(h)) ) (x,x'). 



(14.15) 

This formula proves both statements in the Theorem. 



□ 
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When in Definition (30 1 we introduced the notion of measure change, convergence and smooth- 
ness in the continuum limit of the measure change function G y (y'; t) was not emphasized. How- 
ever, this concept is important, as it is stressed by Girsanov's theorem below, of which we give 
two independent proofs. 

Theorem 23. (Girsanov.) Consider two lattice diffusions on the simplicial sequence A m with 
generators 

(14.16) C m =n{x,t)V h + a{ ^A h , C h = fi(x,t)V h +^ X ^A h , 

respectively. Here fi(x,i), p,(x,t),a(x,t) are assumed to be smooth functions. Then the Markov 
generators C m and C m are related by the smooth measure change with function 

(14.17) 

V cr{x,t) 2 
and the Radon-Nykodym derivative is given by 

(14.18) ^.)=expj| dy t ^—^ dt + 0(h) 

First Proof. Let v(x, t) = p>(x, t) — /j,(x, t) and let a(x, t) — ^t) 2 1 Consider the measure change 
with 

(14.19) G%(x') = e <*,*)&-*) 4 

The off-diagonal elements of the transformed Markov generator are 
(14.20) 

L,[x,x±n,t)- y 2 h 2 2h J ~ 2h 2 2h + 4a(x,t) 2 + 2a(x,t) 2 + [ h 



The diagonal elements instead change as follows 

v(x,t) 2 /j,(x,t)v(x,t) 
2a{x 1 t) 2 a(x,t) 

Second Proof. Let us consider two generators differing by the drift term 



(14.21) £(x, x; t) = C(x, x; t) - - ^ ' j 1 2 ' ; + 0(h). 



(14.22) C(t) = ?^-A h +n(x,t)V h , C(t) = a( ^A h + ii(x,t)V h 

and consider the formula 

E[ e fo o,(x s .s)dx s +b(x s ,s)ds S ^ Xt = y ^ Xo = x j = 

e Wv,t)-4>(x,o)) Pexp ^ jt ds ^ £ ( s ) + fe(s) _ | CT ( S ) V(s) + ( ft j j ( Xj y y 

(14.23) 

To derive Girsanov's theorem, let us ask for what choice of the functions a[x, t), b(x, t) the 

right-hand side equals e tc (x, y; t) up to terms of order 0(h). Since 

(14.24) 

e*(c + b-l-o- 2 a , )e-' i > = ^-A h + (a(x)a(x) 2 + n(x))\7 h +^-a(x) 2 + n(x)a(x) + b(x) + 0(h) 



2 ) 2 v v ' v ' ^ v " 2 

we see that the right hand side equals C up to terms of order O(h) if 
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15. Markov Bridges 



A first simple application of this general framework leads to an extension of the results in the 
previous section to the case of a general Markov generator. 

Consider a Markov process on the simplicial sequence A m — h m Z n [—L,L] C K. According 
to theorem (fli]) and assuming that hypothesis MG1, MG2, MG3, MG4 hold, the symbol of a 
Markov generator in canonical form can be expressed as follows: 



fi, s -/ , sinp/i 2 cosph — 1 
£ m (x,p;t) = i[i(x;t)—- \-a(x;t) — 



Although not necessary for the validity of the calculation and the convergence in the limit h — > 0, 
we restrict to the special case where the generator has the form 

(15.1) C m (x, y; t) = ifi(x; t)V h + a(x; t) 2 A h + h m X m (x, y; t) 

where the functions /j.(x] t) and a(x;t) are smooth in both arguments and X m (x,y;t) > is 
smooth and non-negative for x ^ y and we have that 

(15.2) X m (x,x;t) = - £ X m (x,y;t). 



y=£x£A„ 



Consider a stochastic integral of the form 



It = a(x s ,s)dx s + b(x s ,s)ds 



(15.3) 

where a(x,t) and b(x,t) are smooth functions in both arguments. Let <f>(x,t) — f Q a(y,t)dy. 
The following result holds: 

Theorem 24. (Markov Bridges.) The joint probability distribution function for I t and the 

underlying process on the bridge leading from x to x' is given by 

Pexp y /J L m (s)ds\ (x, J; x', I') 
= I le^-^-^+^lPap ( /' (C m (s) - ipb(s) + l ia 2 a'{s) ~ j>(s) + « m , p (s) + 0(h))ds)(x,x'). 



(15.4) 

where K m p (t) is the operator of matrix elements 



(15.5) n mtP (x, x'; t) = h m X m (x, x'; t) exp f — ip[a(x, t)(x' — x) + <fi(x, t) ~ <fi(x' , t)\ 
f. Consider again a discretization of th< 
C m (x, n; x' , n') — C m {x; x')5 \n' — n — 



Proof. Consider again a discretization of the process It = (Al)nt and the lifted generator 

a(x)(x' — x) 



(AI) 



+ \0n',n+l — On'n) 



(15.6) 

where [a] stands for the nearest integer to a. The partial Fourier transform in the n variable of 
this kernel is 



C mtP (x,x') = ££(:£, 0;x',n)t 



-i(AI)pn 



(15.7) 



= £ m (x; x') exp ( - i(AI)p 



a{x){x' — x) 



(AI) 



(AI) 



( e ^(A/) _ y 



In the limit as (AI) — * we find 
a(xf 



(15.8) 



} A h + (jti(x) - ipa(x) 2 a(x)) V h ~ ipa(x)fj,(x) - -p 2 a (x) 2 a(x) 2 - ipb(x) 

+h m X m (x, x'; t )e- ipa ^'-^ + 0(h) 
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Introducing the function cj>(x,t) defined above, we find that 

e ip ^£ m (x, x'; t) e -^ x '^ = C m>p (x, x'; t) + ipb(x, t)5 xx , - |a(z, t) 2 a'(x, t)S xx , 

(15.9) +h m X m (x , x 1 ; t) exp y — ip[a(x,t)(x' — x) + 4>(x,t) — (f>(x',t)]^ + 0(h m ). 

Hence 

exp (tL m , p )(x, x'; t) 

=ip ^,t)-^',m Pexp ( /„* (C m (s) - ipb{s) + f a(s) 2 a'(s) - j>(s) + K m , p (s) + 0(y)i«J (x,x'). 
(15.10) 



The joint kernel is thus given by (16.14). □ 

16. Abelian Processes in Continuous Time 

This section is based on work in collaboration with Harry Lo and Alex Mijatovic, see (Albanese 
et al. 2006a). 

Abelian processes are path-dependent processes which provide an extension of the notion of 
stochastic integral for which one can extend the Feynman-Kac theorem and Ito's formula in 
Section ([Hi. 

Let A m , m > m be a simplicial sequence, consider a time interval [T, T'\ and a stochastic pro- 
cess described by the Markov generator L m (yi; y 2 ; t). To numerically exponentiate, it is crucial 
in most application examples to first reduce dimensionality by means of block diagonalisations, 
or else the multiplication of the lifted generators would not be feasible. It turns out that this is 
possible for the payoffs of the type above as these two fall in the general class outlined by the 
following definition: 

Definition 43. (Markov Generator.) Let £(j/i; J/a; t) be a Markov generator and let us con- 
sider a lifting of the form 

(16.1) L(yi,ki;y 2 ,k 2 ;t) = L(y 1 ;y 2 ;t)S kl k 2 + A(yi,k 1 ;y 2 ,k 2 ;t). 

where k — 0...K. This lifting is called Abelian if the operators A(yi;y 2 ;t) defined for each pair 
2/i 1 2/2 G A and all times t S [T, T'] as the linear operators of matrix elements 

(16-2) A{yi,y 2 ;t) klM = A(yi, fa; y 2 , k 2 ; t) 

are mutually commuting in the sense that 

(16.3) A(yi,y 2 )A{y 3 , y 4 ) = A(y 3 , yi )A( yi ,y 2 ) 
for all 2/1)2/2)2/3)2/4 £ A an d a ^ t G [T,T']. 

Lemma 4. If all the matrices A(yi;y 2 ;t)k lt k 2 > 2/1)2/2 £ A, ! £ [T, T'] are mutually commuting 
and if furthermore they are all diagonalizable, then they can be diagonalized simultaneously at 
any given point in time, i.e. there exists a time dependent matrix V(k,i;t),i = 0,..K such 
that V{t)~ 1 A{y\, y 2 ; t)V(t) = A(yi;y 2 ;t), where A(yi;y 2 ;t) is a diagonal matrix of the form 
(Hyi;y2;t)8 iu i 2 ) for any t G [T,T'\. 

Proof. Fix a t G [T, T']. If 2/1,2/2,2/3, 2/4 G A and A(yi;y 2 ; t)ip(t) = X(yi; y 2 ;t)ip(t) for some vector 
4>(t), then 

(16.4) A(y 3 ;y 4 ;t)A( yi ;y 2 ;t)ij = A( yi ;y 2 ;t)A{y 3 ;y A ;t)<ip{t) = A( y i;y 2 ;t)A(y s ; yi ;t)ip(t). 

Hence if ip(t) is an eigenvector of A(yi; y 2 ; t) also A(y 3 ; y 4 ; t)ip(t) is an eigenvector of A(yi; y 2 ; t). 
If ip(t) is an eigenvector of multiplicity one, this shows that ip(t) is also an eigenvector of 
A(y3', 2/4; t) for all J/3, j/4. Otherwise, we conclude that the eigenspace of A(y\; y 2 ; t) of eigenvalue 
A is invariant under A(y 3 ,y4) for any j/3,y 4 . Iterating the argument above to this eigenspace 
one can constructively obtain a common set of invariant eigenspaces shared by all the operators 
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■M,Vi'i Vi \ t) for j/i, j/2 € A. Hence, for any given t £ [T, T'], these operators can be diagonalized si- 
multaneously by the matrix whose column vectors span bases for all the common eigenspaces. □ 

Let V(k,i;t),i — 0, ..K be a matrix which diagonalizes simultaneously all members of the 
family of operators T(ki, k 2 \ Vii), so that V(t)~ 1 T(yi; t)V(t) = A(j/i;i) is diagonal. Consider 
the lifted matrix V(fc, j/i; i, j/2; t) = V(k; i\ t)S yiy2 and the transformed lifted Markov generator 

(16.5) (V(t)~ 1 LV(i))(j/i, i\\ j/2, i2',t) = L(j/i; j/ 2 ; i)<5j li2 + A(y 1 ;t) il 

Since this matrix is block-diagonal, the problem of exponentiating it is reduced to the problem of 
exponentiating K blocks separately. This reduces the overall numerical complexity and makes 
it comparable with the complexity of evaluating propagators on K different time points. As 
a further simplification, to exponentiate this block-diagonal matrix it is necessary to hold in 
memory only the blocks, and not the entire matrix. Notice that the blocks are in general 
complex matrices. Hence to fast exponentiate them one needs to invoke the complex matrix- 
matrix multiplication routine sgemm or cgemm of Lcvcl-3 BLAS, depending on whether the block- 
diagonalizing transformation is real or complex valued. 

Example 7. (Stochastic Integrals) Stochastic integrals over diffusion processes provide ex- 
amples of Abelian process. This was already used in the previous section but it is worthwhile 
here to stress the link between the computability of the characteristic functions for the path 
dependent process on a bridge with the block-diagonalizability under partial Fourier transforms 
of the lifted generator. Working directly with Markov generators, one can also generalize the 
results in the previous sections. 
Consider the integral 

(16.6) I t = I(y. , t) = f ( </>(y s ; s) + X (j/ S _ , y s ; a) lim ^ Vs ' s) ~ ^ Vs ~ t] S ~ * } ) ds 



JT \ *i° t j 

where </>(y; s), ip{y, s) and xiyi\ IJ2', s) are continuous functions. Consider the problem of finding 
its distribution on a bridge for the underlying Markov process y t . Introducing the discretisation 

(16.7) I t w (AI)m t 

where mt € [0, ..N] and (A/) > 0, the lifted generator for the pair of processes (yt, (Al)mt) is 
Hyi,m 1 ;y 2 ,m 2 ;t) = L(y 1 ;y 2 ;t)8(m 2 - mi - [(A/) _1 x(j/i; y 2 ; t)(ip(y 2 ; t) - ^(j/i; t))]) 

(16.8) 

Theorem 25. (Abelian Bridges.) The joint kernel of I t and the underlying process on the 
bridge leading from j/i to y 2 is given by 

Peti L ^( yi ,I i; y 2 ,I 2 ;s) = J g e ^-^)p e xp ^J* (C m (s) - ip<j>(s) + K m>p (s) + 0(h)f) (x, x'). 
(16.9) 

where K mtP (t) is the operator with matrix elements 

(16.10) n mtP (y 1 ,y 2 ;t) = £ ro (j/i,j/ 2 ;t)fexp ( - ipx(yi, J/2,*)(V , (jta)*) -^(j/i.*))"^ - l\ 

Notice that this theorem extends the result in the previous section on Markov bridges as it 
includes multifactor processes such as processes with stochastic volatility and, more generally, 
regime switching. 

Example 8. (Double Liftings) There are situations that emerge in practice where one has to 



track two integrals over paths, one of form (16.6) and a similar one of form 
(16.11) I' t = I'(y. , t) ee £ (d>'(y s , s) + x '(y s _ , y. t s) Urn 
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In this case one can introduce a similar approximation 
(16.12) I' t « T t = (AI')n t 

where n t £ [0, ..N] and a double lifting for the generator 
L(y 1 ,m 1 ,n 1 ;y 2 ,m 2 ,n 2 ;t) = 



L(j/i;j/2;*)5( m 2 - mi 



(16.13) 



X(yi;y2;t)(ip(y 2 ;t) - ip(y!;t)) 



(AI) 



S\ n 2 — n\ 



x'(yi;y 2 ;t)W(y 2 ;t)-jj'( yi ;t)) 
(A/') 



s s 

+ -^^^(yi-,t)S mi+ l, m2 Sn u n 2 + +J^<P\yi-,t)S ni + l!n2 5 mum2 . 



Theorem 26. (Multifactor Bridges.) The joint kernel of It, I't o,nd the underlying process 
on the bridge leading from y\ to y 2 is given by the double Fourier transform 

e tL ™(yi,h;y 2 ,I 2 ) = J f^Cwo j <^ e i P '(i' 2 -i[) 

Pcxp ^/ * (£ m (s) - ip<j){s) - ip'(f>'(s) + K m ,p iP '(s) + 0(h))j(x,x'). 

(16.14) 

where n m , P ,p' (t) is the operator with matrix elements 
(16.15) 

K m ,p, P '{yi]y 2 ]t) = L^y^y^t)^- 1 ^ 1 '^^ 2 ^^ 

Example 9. (The Max Process) Consider the path dependent quantity given by the maxi- 
mum attained by a given Markov process, i.e. 

(16.16) K t = K(y t ,t) =ma,x X (y s ;s)ds 

Let's introduce the approximation 

(16.17) K t nK t = ak u 

where a is a constant and k t is an integer value process m t whose paths take values in k t € 
{0...M}. The dynamics for the joint process (yt,k t ) is defined by the lifted generator L on A 
such that: 

(16-18) L(yi, fci; y 2 , k 2 ; t) = L(y 1 ; y 2 ;t)8 klM + A{yi) kl M 5 viV2 

where 



(16.19) A( yi ) kuk2 = 



A if x(yi) > a*i and k 2 = [x{yi)/a] 
otherwise. 



where [a] stands for the nearest integer to a and A > is a fixed number. Typically, A is chosen 
to be a large number as the approximation converges in the limit at A — > oo and a — > 0. A 
direct calculation shows that the operators A(yi) kltk2 commute and hence the maximum of the 
underlying process is itself an Abelian path dependent process. 

17. Abelian Processes in Discrete Time and non-Resonance Conditions 

This section is based on work in collaboration with Manlio Trovato (Albanese and Trovato 
2005) and Paul Jones, see (Albanese and Jones 2007). 

An important class of path-dependent options requires computing the joint distribution of 
the underlying lattice process and of a discrete sum of the following form: 

n 

(i7.i) J t = J(y t ,t) = J2^(yti-^yu;u) 
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where n is an integer, ti = iAT and t = nAT. Suppose that the Markov generator L is time 
homogenous in the interval [to^n] and consider the elementary propagator 

(17.2) Ufa,y 2 ) = Pexp (J ^ L( S )d S ^j (y u y 2 ). 

To find the joint transition probability, one can again discretize the variable Jt so that 

(17.3) J t w (AJ)n t . 

As opposed to lifting the generator as in the continuous time case discussed in the previous 
section, here we lift the elementary propagator itself and form the joint propagator 

(17.4) Ufa,k 1 ;y 2 ,k 2 ) = Ufa,y 2 )S(k 1 -k 2 + [i>fa,y 2 )(AJ)- 1 ]). 

This operator can be block-diagonalized by means of a partial Fourier transform (Albanese and 
Trovato 2005). 

Theorem 27. Consider the Fourier transform operator Ui tP of matrix elements 

(17.5) U it fa;y 2 ) = Ufa, y 2 ) e -^^^). 
Then we have that 

(17.6) (p J] u) (y 1: J 1 ;y 2 , J 2 ) = f (pfj (2/1,2/2) 

This case also falls under a more general class of Abelian liftings. 

Definition 44. (Propagator Liftings.) Let Ui(yi,y 2 ) be a family of Markov propagators 
indexed by i = 0, 1, ...n— 1 defined over the time intervals [tj, tj+i]. Consider a propagator lifting 
of the form 

(17.7) Ufa,ki;y 2 ,k 2 ) = Ufa,y 2 )Qfa,y 2 )k 1 k 2 - 

where k\,k 2 = 0...K . This lifting is called Abelian if the operators Qi(yi, y 2 ) satisfy the following 
commutation condition: 

(17.8) Qfa,y2)Q i {y' 1 ,y 2 )-Qi{y' 1 ,y 2 )Qfa,y 2 ) = 

for all 2/1,2/2,2/3,2/4 € A m and all time intervals i = 0, ...n — 1. 

Let V(k,j),j = 0,..K be a matrix that simultaneously diagonalizes all operators of the 
family Q(yi,y 2 )k 1 k 2 so that V^ 1 Q(yi,y 2 )V = A(yi,y 2 ) is diagonal. Consider the lifted matrix 
V(fc, 2/1; j,y 2 ) — V{k,j)5 yiy2 and the transformed lifted Markov generator 

(17.9) (V-'UiVfaJ^y^h) = Ufa;y 2 )Afa,y 2 )S jlh . 
Since this matrix is block-diagonal, the evaluation of the time-ordered product 

n—l / n—1 v 

(17.10) p ii u = Yn( ri ^ v ) ri 

i=0 ^ i=0 ' 

involves only multiplying blocks whose dimension equals the size of the lattice A m . The reader 
will recognize that the formula in (17.101 is yet another generalization and extension of the 
Cameron-Feynman-Girsanov-Ito-Kac-Martin formulas discussed above. 

The non-singular linear transformation that accomplished the block-diagonalization can be a 
Fourier transform or a more general transformation. Several possibilities are open and the opti- 
mal choice depends on the objective. When using transforms more general than the Fourier trans- 
form, one has to keep present that the simultaneous diagonalization of the matrices Q(yi, 2/2)^1 * a 
above can possibly result in a numerically ill-conditioned problem. An example of this phenom- 
enon arises when one attempts to use a non-homogeneous discretization of the path process 
coordinate. To seize the benefit of inhomogeneous lattices, one needs to be careful when dis- 
cretizing and avoid resonances. 
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Consider a lattice for the path dependent process defined as follows f2 = {o>i, &i, ....ljk } where 
a>i < o-> 2 < ■•■ < W K- Consider the shift operator 7?. of matrix elements 

(17.11) R kl k 2 = (1 -Pk 1 )Sk 1 k 2 +Pfci4 2 ,fci+i, 

if fci < &2 and -Rfcx = 0. Here we assume that Q < pk < 1, for k = 1, ..n. The eigenvalues of 1Z 
are given by the diagonal elements pk = 1 — p k - Diagonalizing this sort of matrix is potentially 
a seriously ill-conditioned especially if the matrix is large and the eigenvalues are degenerate of 
nearly degenerate. For practical application, one thus needs to achieve non resonance conditions 
and keep the eigenvalues pi as far apart from each other as possible. 

Given a non-resonant choice of transition probabilities such as the one above, one can then 
fix a ujq and a Auj > and determine the lattice £1 in such a way that 

K 

(17.12) R k!k 2 ^k 2 = uj kl + Auj 

k 2 =\ 

for all h = l,....K - 1. 

To obtain a non-resonant spectrum, one may choose the lattice u>k so that 

(17.13) io k = u ■ Z k 

where uj$ > and Z > 1. For this discretization to be acceptable, the resulting transition 
probabilities 

u Z k (Z - 1) 

ought to be between and 1. Hence Z must satisfy the additional constraint 
(17.15) Z(Z -!)>—. 

LOo 

We interpolate between the transition kernels 1Z to obtain the relevant kernels of all possible 
values of t/j(xi, X2). We construct the relevant kernels as follows: 

Q(vi,V2) = Myi,V2)] + 1 - </%i,2/ 2 ))^ ( ^ 2)1 + Mift.ift) - [^{yxM])^ mvuV2)]+ \ 

(17.16) 

for all yi and t/2 in A m , where [ip{yi,y2)] represents the integer part of the functional ip at 
(2/1 1 2/2)- In the representation in which the operator 1Z is diagonal, the operators Q (2/1,2/2) are 
also diagonal and the lifted propagator is block-diagonal. 

18. Univariate Moment Expansions on Bridges 

This section is based on work in collaboration with Adel Osseiran, see (Albanese and Osseiran 
2007). 

A second methodology that applies to most Abelian path dependent options is is based on 
obtaining only a few moments of an Abelian process on any given bridge with respect to the 
underlying Markov process, as opposed to reconstructing the entire conditional distribution. 

Consider a time interval [T, T'\ and a Markov generator L(j/i, 7/2; t). Consider again the 



integral it in (16.6 1. Let's introduce the following one parameter family of deformed Markov 



operators parameterized by e e R 

(18-1) Myi, 2/25 £) = F(yi,y 2 ;i) + e-V(y 1 ,y 2 ;t) 

where 

(18.2) V(y!,y 2 ;t) = <f>(yi;t)5 yit y 2 + L(y 1 ,y 2 ;t)x(yi,y2;t)(^(y2;t) - i>{yi\t)) 
Theorem 28. (Dyson Formula.) We have that 

(18.3) 



de 



=0 



Pexp (J £ e («)<ifl)(»i,lft) = E T [r t l 5(y t - y 2 )\y T = 2/1] 
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Proof. Consider Neper's formula for the propagator 

(18.4) Pe^(^j\ e (s)ds^ = Jim^Pn( 1+ ^r"( L ^) + *))) 

where ti = T + j^. By collecting similar powers in e, one finds Dyson's formula 
(18.5) 



Pexp ^ J h e (s)dsj = exp((i - T)L+ 

(18.6) e J ds 1 (^ (si - T '> L V(s 1 )e (t - Sl ^ + 

(18.7) e " /* ds i- f g ds ™ ^ Sl - T '> L V(s 1 )e^- Sl '> L ....V(s n )e^- s ^ L ^ . 
The time-ordered integrals above are proportional to conditional moments, i.e. 

/ft \ oo £ n 

(18.8) Pexp ^ L e (s)dsj ( yi ,y 2 ) = £ -£ T [/ t "% - y 2 )|2/T = . 

Here, the factorials originate from the time ordering. □ 
A technique which is numerically stable in many situations of practical relevance is to nu- 
merically differentiate with respect to e the deformed propagators Pexp ^ L e (s)ds^j (j/i, j/2) 

and evaluate at e = 0. This technique can be used to obtain the first moments of It on any 
given bridge for the underlying Markov process. In most applications, we find that two moments 
suffice to extrapolate the probability distribution function to sufficient accuracy. To do so, it 
is convenient to choose from among the probability distribution functions which are analyti- 
cally tractable. For instance, starting from the first two moments only, one can use either the 
log-normal or the chi-square distribution. 

The standard chi-square distribution is given by 

where a is the number of degrees of freedom. The first and second (raw) moments of this 
distribution are 

E[x] = a, E[x 2 ] = a(a + 2) 

To match the pre-assigned first and second moments mi, 777,2, respectively, one can pass to the 
new variable 

7771 

£ = — x 

a 

and chose 

2m\ 



a = 



m 2 — ra\ 



Let £ be a log-normally distributed random variable with probability distribution function 



The first two moments are given by 

E[x] = e»+° 2 ' 2 and E[x 2 ] = e 2 ^ 

The parameters \i and a can be reconstructed from two pre-assigned first and second moments, 
mi and 7772 respectively, by setting 

[i = log ( J_ ) and er 2 = log ' 



/m 2 
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The log-normal and the chi-square distributions allow one to match the first two moments. 
More general distributions available in closed form allow one to match higher moments also. The 
Pearson Type III distribution for instance has a probability distribution function given by 

f(x) - 1 ( x ~ a Y 1 e -(x-o)/6 
nX) -bT(p) { b ) 

and extends over the half line [a, +oo). The special case of this distribution when a = 0, b = 2 
and p is half of an integer, gives the Chi-Squarcd distribution. In general, the moments are given 

by 

E[x] = a + bp 

E[x 2 } = (a + bp) 2 + b 2 p 

E[x 3 ] = (a + bp) 3 + 3b 2 p(a + bp) + 2b 3 p 

and matching these with our pre-assigned moments mi , m 2 and 7713 and computing the values 
of a, b and p in terms of these moments we get 

2(m 2 -ml) 2 

a = mi — = 

to 3 + 2m\ — Jmim 2 

m^ + 2m\ — imim 2 
2(m 2 - m\) 

4(m 2 - m?) 3 

P = 

(m 3 + 2m\ — 3mim 2 ) 2 

Example 10. (Volatility contracts) As an example, consider variance and volatility swaps. 
Variance swaps of maturity T" and time of issuance T have a payoff given by 

min ^jT v(y s _ Q ,y s )ds, f ■ SR 2 ] - SR 2 



where SR is the swap rate and / is a factor, a typical value being / = 6.2. Here, v(yi,y 2 ,t) is 
the instantaneous variance defined as follows: 

(S(V2)\ 

\S(yi)J 

if yi is such that x(yi) ^ 0. Otherwise, if x(y\) = then v(y\,t) = 00. The variance swap is 
said to be at equilibrium if its price is 0. The payoff of a volatility swap is 



(18.9) v(y u y 2 ) = jC(yi,y 2 ;t)log 2 



min I v(y s ,s)ds,SR ] — SR 

Here, SR is the volatility swap rate. 

To price these contracts it suffices to evaluate the distribution of realized variance on a bridge, 
i.e. of the functional 

i-T' 

Is 

It 

By approximating the distribution of RV(y) with the chi-squared distribution, and using the 
CDF 

7(f,f) 



T' 

(18.10) RV(y 2 ) = S(y T '-y2) J v(y s ,s)ds 



F(x;a) = 



where T(z) and 7(^,0) are the gamma and incomplete gamma functions respectively, i.e: 

poo pa 

T(z)= s'-^-'ds, 7(3,0)=/ s'^e-'ds 



we find 

£ t [min(W, W max ) <5(y T , - y 2 )] = ^ [if (1 - F(X; a)) + aF(X; a + 2)] 
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and 



Et 
where 



min(VKV, \J RV max ) 5(y T , - y 2 ) = 



K(1-F(K; o)) + 



r(f) 



2mi(yi,y 2 ) ~ ~, s 0(3/1,3/2) DT/ 

a = o(j/i,J/2) = 7 : 7 rrj, K = K(y 1 ,y 2 ) = ? r-KVW 

m 2 (2/1, 2/2) - mi (3/1, 3/2 r mi (3/1,^/2) 

Since the dependency on the swap rate in both cases is non-linear, an exact calculation requires 
using a root finder. 

Using instead the log-normal distribution to extrapolate, we find 

E t [mm(RV, i?U max ) 5(2/ T , - y 2 )] 

= e^' 2 M ( WVn^)-»-° 2 /^ + i? T/ max TV ( '°^ w r >-" ) . 

(18.11) 

For volatility swaps instead, we find 



mm(\/KV, y/RVjnax) 5(y T > - 3/2) 



(18.12) 

where it and cr are specified above in terms of mi and m 2 . 

Finally, in case the Pearson Type 3 distribution is used, we can express these expectations in 
terms of the chi-square cumulative distribution function F(x\ a) as follows: 



E t [mm(RV, RV max ) 8{y T , - y 2 )] = (a + bp- RV max )F chi [ 2p, 2^^- 



RV ml >x - a\ p e 



(18.13) j . + (2RV max -a-bp). 

19. Bivariate Moment Expansions on Bridges 

This section is based on work in collaboration with Adel Osseiran, see (Albanese and Osseiran 
2007). 

Consider a time interval [T, T'] and a Markov generator L(yi, 3/2; t). We are interested in the 
joint distribution on any given bridge for the underlying process of two stochastic integrals 



(19.1) i ltt = i(y.,t) = J ^Mys-,s) + xi(ys-o,y s -,s)iim Mvs ' s) MV'-*' 8 

and 

(19.2) i 2 , = /,:,../) - jf (Uvs; s) + My s -o,y s ; s) Um " "),/„. 

To handle this problem using the moment method, we introduce the following two-parameter 
family of deformed Markov operators parameterized by e, e' e M 

(19-3) L eii£2 (3/i,3/ 2 ;i) = L(y 1 ,y 2 ;t) + ei Vi (3/1, 3/2 ; i) + e 2 V 2 {yi, 3/2; i) 

where 

(19.4) ^1(3/1,3/2;*) = ^1(2/1; i)5j/ ll!/2 + L(j/i, ?y 2 ; *)xi(2/i, 2/2; t){ipi{y 2 ; t) - i>i{yi;t)) 
and 

(19.5) V 2 (y 1 ,y 2 ;t) = fo{yi;t)5 yity2 + L(y 1 ,y 2 ;t)x2{yi,y2;t)(il> 2 (y 2 ;t) - ^2(2/1! *)) 
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Theorem 29. (Dyson Formula.) We have that 



(19.6) 



Qn 1 +n 2 

de^de? 



Pcxp ( f L ei , e2 (s)ds)( yi ,y 2 ) = E T [l^ t I%5(y t - y 2 )\y t = Vl \. 

Proof. The proof is a simple extension of the proof in the univariate case and will be left to the 
reader. □ 

Consider a bi-variate log-normally distributed variable 

Yi\ f log^i) 
Y 2 J \ log(X 2 ) 

whereby X x and X 2 are correlated normally distributed variables of joint distribution 
(19.7) f{x u x 2 ) 



N 



Mi 

M2 



°~2 



exp ■ 



1 

W^J 2 



log X X - p 1 



-2p 



2-KO\0 2 \]\ - p 2 X\X 2 

(log si - /ii)(loga;2 - p, 2 ) (\ogx 2 -p 2 



o\o 2 



o- 2 



where p is a correlation parameter. Both X\ and X 2 are log-normally distributed with 

E[Xi] = e ^+^/ 2 , and E[Xf] = e 2 ^ +2rT > , i = 1,2. 
Matching these with the pre-assigned moments, and solving for p,- L and cr, we find 



Hi = log \ E 



(i) 



(4->y 



and 



log £ 



(i) 



Moreover, the mixed term is 

£ [XiX 2 ] = E [e Yl+Y2 ] = e ^i+M2+i(^+<T 2 2 +2^i<r 2 ) = E [ Xl ]E[X 2 ] e" CTlC 
which gives p (in terms of the pre-assigned moments): 



(19.8) 



P = 



<7 1 a 2 



log 



E 


'r(l) r (2)' 




E 






r(2) 



Example 11. (Conditional Variance Swaps) The payoff of a conditional variance swap is 
given by the ratio 



(19.9) 



$t v(y s ,s)l(L<S(y s ) < H)ds 

™ ail 



£ 1(L < S(y a ) < H)ds 



To apply the moment method to this case, the first thing to note is that essentially we are 
modelling the two integrals appearing in the payoff at the same time. We are going to need a 
Bi-variate distribution to do this. Firstly let's write: 



It 1] =J T v(y s ,s)l(L<S(y s ) < H)ds 



r(2) 



/; 



1(L < S(y s ) < H)ds 



and in order to compute the expectation 
(19.10) E 
we'll need the following expectations: 



r(2) 



E 



r(i) 



E 



r(2) 



E 



, and E 



r(l) r(2) 
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To tackle this problem consider the operator 
(19.11) 
where 



(19.12) 

and 
(19.13) 

We have that 
and 



C eu€2 (yi,y 2 ) = C{y!,y 2 ) +e 1 (j>(y 1 )S yi y 2 + e 2 ?/%i)<W 
0(y ls t) = ^£( yi ,y 2 ;i)log 2 (^\) 1(L < S( Vl ) < H) 



d_ 
del 



ei=0 



8e\ 



\S(yi)J 
i>(y 1 )=l(L<S(y 1 )<H). 
eW-W*i,-»\y 1 ,y 2 )=E 

e^ t '- t ^—\y 1 ,y 2 )=E 



r(i) 



Vt = ViiVt' = 2/2 



£1=0 



similarly (but with respect to e 2 ) for 



E 



n 



(2) 



yt = yuvv = v-2 



and E 



yt = yi,y t = 2/2 



yt = yuw = V2 



The joint expectation will involve the mixed derivative: 

d 2 



(19.14) 



E 



r(l) r(2) 
1 t 2 t 



de\de 2 



((t'-t)£ ei>e2 ) 



(2/1,2/2) 



and once computed, we make use of these expectations to match a bivariate distribution. For 
simplicity of notation we leave out the conditional part of these expectations, noting that all the 
moments we obtain are conditional to the initial and final points. 

let us notice that 



To evaluate the expectation E 

E 



r(l) / r(2) 



Xn 



= E [e yi - F2 ] = E 



where Y 2 = — Y 2 is also be normally distributed ~ N(—fjL 2 J cr 2 )- Hence 

~X 1 ' 



E 



X, 



; Ml + 5 <T l g-/*2 + 3°2 e -pfflCT2 



and the expectation (19.101 is given by 

E 



(19.15) E 



(i) 



r(2) 



E 



(i) 



E 



(2) 



E 




E 


7(2)" 


E 


V(l)r(2)" 
1 t 1 t 





E 



(i) 



E 



E 



r(2) 



E 



r(l) r(2) 
1 t 1 t 



moreover 



so 



(19.16) 



E 



x 2 



E 



^2n 1 +2o 1 e -2fj, 2 +2cr 2 e -4pa 1 cr 2 _ 
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This double integral will need to be evaluated numerically. As the first one has finite bounds and 
the second has an infinite upper bound it would make sense to use two types of Gaussian quadra- 
tures: a Gauss-Legendre quadrature on the inner integral and a Gauss-Laguerre quadrature on 
the outer one. 



20. Block Factorizations 

Although most path dependent processes emerging in applications to Finance are Abelian, 
some aren't. If the Abelian property fails, the propagator cannot be block-diagonalized by the 
methods discussed above and also moment methods generally fail. In several relevant situations 
one can still achieve a numerically viable framework by using block-factorizations instead of 
diagonalizations. 

Let Ui(yi, y 2 ) be a family of Markov propagator indexed by i = 0, 1, ...n — 1 defined over the 
time intervals [tj,tj + i]. Consider a propagator lifting of the form 

(20.1) U i (y 1 ,k 1 ;y 2 ,k 2 ) = U i (y 1 ,y 2 )Q l (y 1 ,y 2 ) klk2 . 

where ki,k 2 = 0...K. 



Definition 45. (Block- Factorization.) We say that the propagators in (20.1) admit a block- 
factorization if they can be represented in the form: 

(20.2) Ui = IU ■ (Ui <8> I) 

Here TLi,i = 0, ...n — 1 is a family of permutation matrices, i.e. matrices with the property that 
for each pair (t/i, k\) we have that 

(20.3) ll t (y 1 ,k 1 ;y 2 ,k 2 )=0 

for all pairs (y 2 ,k 2 ) except for one single pair (Y{yi,ki),K{yi,ki)) for which we have instead 

(20.4) ILi( yi , fci;y( tfl , h),K( yi , h)) = 1. 



Furthermore, the tensor product notation in equation {20.2) stands for the operator with the 
following matrix elements: 

(20.5) (U i ®l)(y 1 ,k 1 ;y 2 ,k 2 ) = Ui( yi ,y 2 )S klk2 . 

If a block-factorization exists, then an efficient backward induction algorithm can be setup. 
Namely, if v is a vector, then we can value the following path ordered products 

n— 1 n 

(20.6) Vi = P Yl Ujv n = i 3 []n j - (Uj <g> I) v n 

j=i j=i 

iteratively in i from i — n — 1 down to i = 0. At the first step one applies the operator (U n -i <8l) 
to v n , followed by the permutation II n _i to obtain u n _i and then iterate. Due to the tensor 
product structure of the first operator, it is convenient to apply by regarding the vector as 
a matrix of elements Vi(y,k). This representation makes it possible to leverage on numerical 
efficiencies and to re-interpret a high dimensional BLAS level 2 matrix- vector multiplication as a 
lower dimensional BLAS level 3 matrix-matrix multiplication. Applying a permutation operator 
is then quite straightforward from the numerical viewpoint and is an operation whose complexity 
scales linearly with respect to the dimension of the vector v. 

Example 12. (Snowballs) Here we follow (Albanese 2007a). Consider the case of a valuing a 
snowballs for which the structured coupon at time Tj = T + (AT)i has the following form: 

(20.7) C Ti = /C Ti _ 1 +$ i (j/T < _ 1 ) 

where the factor / is a fixed parameter and ^(y^-i) is a given function. Since the coupon 
amount at a given time affects the coupon amou/nt in the next period, the process is not 
Abelian, in fact it is path-dependent. However, block-factorizations are still possible. 
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Consider discretizing the coupon dimension in intervals AC, so that a generic coupon can be 
approximated as follows 

(20.8) C Tz = (AC)fc Ti 

where = 0, 1, ....N — 1 is a discrete time, integer value process. A strategy to implement the 
backward induction scheme is to organize the payoff function v(y, k) at maturity in a matrix 
with TV columns, each one indexed by the state variable y and representing the price conditional 
to the discretized value of the previous coupon paid. One can then iterate backward by applying 
the pricing kernel to this matrix of conditional pricing functions. After each step in the iteration, 
one needs to reshuffle the pricing functions in such a way that the conditioning relation on each 
column is satisfied. More precisely, we set 



(20.9) A) - Yl Ui(yi,V2)vi(y2, [(/(AC) • h + ^( yi )) / (AC)}) 



(2/1, *0 



where 

(20.10) Ui(y u fc i; y 2 , k 2 ) = S(y 2 - Vl )8{k 2 - K{ Vl M)) 
and 

(20.11) K( Vl , fcO - [(/(AC) • h + ^( yi ))/(AC)}. 

Notice that this backward induction scheme can accommodate callability provisions. 

Example 13. (Soft Calls) Soft calls provide another interesting example of a non-Abelian 
process. In this case, a call decision depends on the fraction of time the process spends above a 
given barrier level during a fixed time period prior to the decision itself. By restricting times at 
which one makes observations to the discrete sequence of time points ti — T + iSt,i = 0, ...M, 
the pricing function takes the form 

(20.12) Wj (y,-?) 
where 

(20.13) ^ = (<7 1 ,...,a JV )e{0,l} Jv . 

Here the variable Oj equals 1 if yti_ N+j satisfies the barrier condition and otherwise. We model 
this dependency by means of a generic function notation 

(20.14) aj = Vi{y ti _ N+j )- 

The sequences ~a can be arranged in an ordered sequence of 2 N integers so that the data structure 
Vi(y, a ) appears as a matrix with 2^ columns. A backward induction scheme involves evaluating 

(20.15) u»-i(yi, = X! u i(yi,y2) v i(y2,$i{~v, yi)) 

V2 

where 

(20.16) $iCCT,2/l) = {CT2,-0-AT-l,£(*i,*»-l,2/l)} 

This form is block-factorized as 



(20.17) Wi_i(i/i,"?) 



IT • (Ui ®l)vi 



(2/1, *0 



where TE^ is the permutation operator of matrix 

N-l 

(20.18) IIi(yi, ^i; 2/2, "? 2 ) - S{y 2 - yi)S((^ 2 ) N - S,^)) [] 5(tf 2 )j - 

j=i 

More general types of useful block-factorizations can be imagined, although they might be 
numerically less efficient. An example is set forth by the following definition: 
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Definition 46. (State Dependent Block-Factorization.) We say that the propagators in 
(20.1) admit a state dependent block-factorization if they can be represented in the form: 



K 



(20.19) Ut = it • Ui, 



k=l 

Here IT,z = 0, ...n— 1 is a family of permutation matrices as above. The direct sum notation in 
equation (20.19) stands for the operator with the following matrix elements: 

K 

(20.20) (0 U iM> )(y 1 ,k 1 ;y 2 ,k 2 ) = U t . kl (y 1 ,y 2 )S klk2 . 

ko=l 

A backward induction in this case is still a lower dimensional problem but since the lower 
dimensional propagator is state dependent this reduces to a sequence of matrix-vector multi- 
plications as opposed to a single matrix-matrix multiplication. The scheme is thus numerically 
less efficient. Still it would be useful in cases for instance such as GARCH type models and 
extensions, see (Engle 1982). 



21. Dynamic Conditioning 

This section is based on work in collaboration with Alicia Vidler, see (Albanese and Vidler 
2007). 

Modeling correlations between several processes has sometimes been considered a problem 
affected by the so-called curse of dimensionality which causes the numerical complexity to explode 
exponentially see (Bellman 1957). This motivates the recourse to Monte Carlo methods while 
calibrating using closed form solutions. 

Here we deviate substantially from this framework and aim at building models viable from the 
engineering deployment viewpoint and specified semi-parametrically or even non-paramctrically 
without any restriction imposed by analytical tractability. The method of fast exponentiation 
allows one to calibrate single name models without recurring to analytic solvability thanks to 
providing smooth sensitivities. Although Monte Carlo schemes could in principle be employed 
in a model constructed non-parametrically and calibrated with operator methods, in many cir- 
cumstances important for applications this can be avoided by using the technique of dynamic 
conditioning discussed in this section. This technique can be regarded as a dynamic copula 
whereby the single factor marginal distributions are preserved and the numerical complexity 
grows linearly with the number of factors, similarly to what happens with Monte Carlo meth- 
ods. But, unlike Motecarlo methods, dynamic conditioning is numerically noiseless as there is 
no sampling, no variance reduction scheme is needed and even features such as callability are 
numerically treatable. 

Dynamic conditioning is based on constructing a hierarchy of conditioning relations. At 
the base of the hierarchy we have continuous time lattice models for single factor marginal 
distributions. Next, one introduces a binomial process for each risk factor to condition the 
corresponding Markov chain. Next one finds a hierarchy of binomial processes to condition the 
former conditioning lattices. See figure Q. 

To describe the technique of dynamic conditioning in detail, consider a particular risk factor 
described by a lattice process whose filtered probability space is engendered by a Markov propa- 
gator U(yi, t\\ t/2j ti)- For each such single factor, we introduce a conditioning binomial process 
/it 6 Z, which is constant over the time intervals [Tj, Tj+i) where Tj = T + iAT, i = 0, 1, 2, ...N 
and AT = (T' — T)/N. At initial time we set hx — while for all z > we have that 
— hx i _ 1 = ±1. The elementary propagator for the process h t across neighboring time points 
T-is 



(21.1) 



V(hi,Ti; hi+x,T i+1 ) = qi(h l ,i)S hl +i.h i+1 + q-i(hi 7 i)S h ^ 1 , hz+1 
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Figure 4. Dynamic conditioning scheme: a binomial process conditions each 
risk factor and is in turn conditioned by another binomial process. The latter 
is also conditioned by a hierarchy of binomial processes. 



where qi{hi,i),q-\{hi,i) > and qi(hi,i) + = 1. On a general time interval [Ti,Tj], 

we have that 

3-1 

(21.2) V(h l ,T l ;h J ,T J )= £ J] V(h Tk , T k ; h n+1 , T k+1 ). 

ht-.h^—hijhT^—hj k—i 

Next we define a lifted conditional propagator U (j/i, hi,Tf, yj, hj, Tj) in such a way to preserve 
unconditional transition probabilities from the starting time at T, i.e. so that 

(21.3) ^(w. 0: r ; w. ^-» T i) = ^G/o. r ; ^,7)). 

To do so, one strategy is to define two propagators for each node (h,Ti), namely to choose a 
pair of operators Ux{y il h i ,T i \y l j rl ,hi + l,T i+ i) and t/„i(y;, hi, T t ; y i+1 , hi - l,T i+ i) so that 

U(yi,Ti;y i+1 ,T l+1 ) = q^h^ijU^yi, h,Ti, y i+1 , h + l,T i+1 )+ 

q-i(h,i)U-i(y l ,h,T i ;y l+ i,h - l,T i+1 ). 



The operator U(yo,0,T;yj,hj,Tj) satisfying ( |2L3j is defined as the following sum: 
U(ya,ho,T Q ;y j ,h j ,T j )= ^ II ^-^^(^-i)* 

ht:h T =0,h T .=hj yi,..yj-i k=l..j 

Uh Tk -h Tk l (yk~i,hT k _ 1 , Tfe_i; y kl hT k ,T k ). 

Since the operators U\ and U-i do not commute with each other, if not in trivial situations, 
each path h t in the summation on the right hand side of this equation represents a different oper- 
ator. In many situations one can however avoid the numerical complexities of a non-recombining 
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tree by finding modified versions of such operators so that the kernels 

U h .(yo,yk)= n q hTk ^ hTk _ i (hT h _ 1 )U(y k - 1 ,hT k _ 1 ,T k -i,y k ,h Tk ,T k ) 

Vl,--Vj-l k=i..j 

are all equal to each other, for all paths ht : hx = 0, = hj and at least one single fixed 
starting point yo = yo- This will be referred to as conditional recombination property. 

The reasons why we focus on the initial point yo — yo only are manifold. Firstly, in applica- 
tions one needs to condition marginal distributions only to spot values, as the price of options 
in the hypothetical case asset prices were different is not known. Secondly, if we insisted on 
the same property being valid for all initial starting points we would end up with a seriously 
ill posed problem of difficult solution. Because of these reasons, we settle for the more modest 
objective of conditional recombination. 

Namely, on each node (hi,T{) with i > we define 

(21.4) U(yi-i,hi-i,Ti-i;yi, hi,Ti) = U hi - hi _ l /li-^Tj-i; j/j, ft,,, Tj) 
in case hi = ±i. Otherwise, we determine this operator so that 

U{yo,ho,To\yi,hi,Ti) = 

Y & (y , h , T ; y i _i,hi- l ,T i _ x )U(y i -i,h i -i,T i _ 1 ; y h hi, T 4 ) 
SK-ieA 

for all yi £ A and a fixed yo. This can be achieved in more than one way. As a guideline, it is 
advisable to deviate the least possible from the definition (21.41. 

As a next step in the construction, consider a second binomial process c t which is picccwisc 
constant across neighboring time points Tj . The dynamics of Ct is by construction correlated to 
that of h t . More precisely, the propagator for the joint process is 

W(hi,Ci,Ti;h i+ i,Ci + \,T i+ i) = 

c l++(hi,Ci, i)5hi+i,hi+i^ci+i,c i+l + Q+-(hi, Ci,i)Sh i+ i t h i+1 5 Ci _i tCi+1 

+ 1-+(hi,Ci, i)8 hz -\,h z+1 5 Cz +i,c t+1 + q — (hi, c i ,i)8 hi -i,h i+1 5c i -i,c i+1 

where q±±(hi, Cj, i) > and 

(21.5) q ++ (hi,Ci,i) + q + -(hi,Ci,i) + q- + (hi,Ci,i) + q — (hi,Ci,i) = 1. 

Due to the conditional recombination property, if yo — yoi then the conditional propagators 
resum with a simple formula 

Ci+i,Ti + i)x 

hf.h T =0,h T j=hj yi,..yj-i k=l..j 

qh Tk -h Tk _ 1 (hT k - 1 )U hTk -h Tk _ l (yk-\,hr k _ x ,T k -i; y k , h Tk ,T k ) 
= w (°> 0. T o; hj , Cj , Tj)U(y , 0, T ; Vj , hj , Tj) 

hj 

= U(yo,Q,To;yj,Cj,Tj). 

As a next step, the construction above is repeated for a number of different risk factors 
associated to N lattice processes y[ a \ where a — 1, ..N all correlated to the process ct, possibly 
in different ways. Then, conditioned to starting all processes at fixed lattice points y^ = y^ 
and conditioned to ctj — Cj , the multi-factor propagator factorizes into the product of conditional 
single factor propagators 

(21-6) J] ^o a) .0,T ; tf j a) ) c, ) T J ) 

a=l...N 

This is the key formula which we use to correlate processes while ensuring that numerical com- 
plexity increases only linearly with the number of factors N. 
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The construction can obviously be iterated and we can have a whole hierarchy of binomial 
processes conditioning each other to model multi-factor correlations. 

22. Conclusion 

We reviewed a new framework for Mathematical Finance and the theory of stochastic pro- 
cesses based on operator methods. The framework grew over the years through applied research 
and solving specific concrete problems in derivative pricing. As wrong ideas were weeded out, 
a coherent framework emerged and is now summarized in this review paper. From the nu- 
merical viewpoint, the methods rely on the ability to multiply efficiently large matrices as the 
main computational engine. Technologies for matrix manipulation are currently being devel- 
oped at great pace. The emerging multi-core GPU chipsets will soon provide affordable teraflop 
performance for algorithms based on matrix-matrix manipulations such as ours. To both set 
the theory of stochastic processes on firm theoretical grounds and justify the empirically ob- 
served convergence rates, we derive pointwise kernel convergence estimates of a new type which 
explain the observed smoothing properties of the method. We introduce the notion of Abelian 
process to generalize the concept of stochastic integral and extend the classic Cameron- Feynman- 
Girsanov-Ito-Kac-Martin theorem in multiple directions. This theorem is here reinterpreted as a 
block-diagonalization scheme for large matrices corresponding to Abelian processes. We outline 
solution methods for Abelian path dependent options, a class we identified and which comprises 
the great majority of the path-dependent options currently traded. Important non-Abelian 
processes are also considered and discussed by means of a weaker but still effective method of 
block-factorizations. Furthermore, we also discuss a method for dynamic conditioning which 
applies to correlation products such as baskets and hybrid derivatives. 
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